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ABSTRACT 


A  local  geoid  model  to  predict  the  geoid  heights  in  the  vicinity  of  Monterey 
Bay,  California,  was  developed  to  use  Global  Positioning  System  (GPS)  differential 
positions  and  known  Mean  Sea  Level  (MSL)  with  the  method  of  collocation.  The 
local  geoid  models  were  based  on  Rapp's  360  degree  x  360  order  global  geoid 
model  determined  from  gravity  measurements.  Control  data  were  adjusted  by  least 
squares  to  solve  for  the  parameters  in  the  local  geoid  model.  Also  studied  were 
factors  that  affected  the  GPS-measured  ellipsoid  height  differences.  These  included 
(1)  comparing  GPS  differencing  solutions,  (2)  standard  error  of  GPS  observations, 
(3)  corrections  for  surface  meteorological  values,  and  (4)  observation  durations  for 
GPS. 

The  data  used  in  this  research  were  taken  from  GPS  measurements  on  the 
campus  of  the  Naval  Postgraduate  School  (NPS),  an  area  about  100  m  x  630  m  and 
in  an  area  approximately  15  km  x  33  km  near  Monterey,  California.  The  time 
period  was  from  February  5,  1988,  to  May  12, 1988. 

The  accuracy  of  the  predicted  geoid  heights  is  ±  2  cm  if  a  six-parameter  model 
is  used  for  the  large  area,  and  ±  2  to  10  mm  if  a  five-parameter  model  is  used  for  the 
NPS  campus. 
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The  Global  Positioning  System  (GPS)  is  able  to  establish  precise  relative 
positions  in  the  World  Geodetic  System  of  1984  (WGS  84).  A  Trimble  GPS 
receiver,  which  has  the  capability  of  measuring  carrier  phase,  was  used  to 
determine  the  vector  base  line  in  space.  The  components  of  the  base  line  are 
expressed  in  terms  of  cartesian  coordinate  differences  (AX,  AY,  AZ)  [Remondi, 
1984].  These  vector  base  lines  can  be  converted  to  distances,  azimuths  and  the 
ellipsoid  height  differences,  Ah,  relative  to  the  WGS  84  Ellipsoid. 

The  results  of  several  tests  and  operations  have  clearly  shown  that  GPS 
survey  methods  can  replace  conventional  horizontal  survey  methods. 
Comparable  accuracies  have  also  been  achieved  for  GPS-derived  ellipsoid  height 
differences.  Ah.  The  problem  of  converting  these  ellipsoid  height  differences. 
Ah,  to  orthometric  height  differences,  AH,  remains  to  be  resolved.  For  example, 
in  engineering  surveying  applications  WGS  84  coordinates  must  be  transformed 
to  the  North  American  Datum  of  1983  (NAD  83)  system.  The  GPS  obtains 
ellipsoid  heights,  rather  than  the  orthometric  heights;  the  geoid  height,  N,  must 
be  calculated  to  obtain  the  latter.  One  of  the  problems  in  this  transformation  is 
the  accurate  determination  of  the  local  geoid  height,  N. 

For  GPS  survey  applications,  a  geoid  model  should  provide  geoid  heights 
with  an  accuracy  commensurate  with  that  of  the  ellipsoid  height,  h,  so  the 
accuracy  of  the  derived  orthometric  heights  is  not  reduced.  In  the  future 
differential  positioning  will  be  widely  used  in  GPS  surveying,  and,  therefore, 
only  the  geoid  height  differences  between  stations  will  be  required. 

Geoid  height  computation  techniques  include  the  following: 

(1)  Geoid  height  differences  in  the  U.S.  can  be  determined  from  gravity 
data  and  the  Stokes'  integral  method,  or  from  astro  gravimetric  data  and  least 
squares  collocation  methods.  These  methods  lead  to  uncertainties  that  are 
typically  1  to  10  cm  for  distances  of  less  than  20  km  and  5  to  20  cm  for  distances 
between  20  to  50  km  [Zilkoski,  1988] 

(2)  A  comparison  of  a  data  set  from  GPS  with  gravimetrically  determined 
geoid  heights  using  least  squares  collocation  techniques  shows  discrepancies 
between  the  two  data  sets  of  about  ±  2  cm  for  a  maximum  intersection  distance  of 
approximately  50  km  [Denker  and  Wenzel,  1987]. 


K 
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(3)  Mean  gravity  anomalies,  deflections  of  the  vertical  and  a  geopotential 
model  calculated  to  degree  and  order  180  have  been  used  to  determine  geoid 
heights  in  the  area  bounded  by  (34°  <  <p  <  42°,  18°  <  X  <  28°)  [Tziavos,  1987]. 
The  method  used  was  that  of  least  squares  collocation.  By  using  empirical 
covariance  functions  for  the  data,  suitable  weighting  functions  for  the  different 
sources  of  observations,  and  the  optimum  cap  radius  around  each  point  of 
elevation,  an  accuracy  better  than  ±  0.60  m  was  obtained  for  geoid  heights. 

The  main  objective  of  this  thesis  is  to  use  a  model  that  predicts  N  in  a  region 
near  Monterey,  California,  from  the  GPS  differential  positions  and  known  mean 
sea  level  (MSL)  using  the  method  of  collocation.  Also  studied  were  factors  that 
affected  the  ellipsoid  height  differences,  Ah,  obtained  from  GPS  measurements, 
such  as  comparing  GPS  differencing  solutions,  standard  error  of  GPS 
observations,  corrections  for  surface  meteorological  value  ,  and  observation 
durations  for  GPS  measurements.  The  data  used  in  this  research  were  taken 
from  GPS  measurements  on  the  campus  of  the  Naval  Postgraduate  School  (NPS) 
an  area  about  100  m  x  630  m  and  in  an  area  approximately  15  km  x  33  km  near 
Monterey,  California.  The  time  period  was  from  February  5,  1988,  to  May  12, 
1988. 

The  accuracy  of  the  predicted  geoid  height  is  ±  2  cm  if  a  six-parameter 
model  is  used  for  the  large  area,  and  ±  2  to  10  mm  if  a  five-parameter  model  is 
used  for  the  NPS  campus. 
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n.  GEOID  HEIGHT 


A.  THE  GEOID 

Surveyors  and  engineers,  in  most  cases,  are  interested  in  the  orthometric 
height,  H,  as  measured  above  the  reference  surface  of  the  geoid.  The  ocean  is 
considered  to  be  freely  moving,  homogeneous  and  only  subject  to  the  force  of 
gravity.  When  a  state  of  equilibrium  is  achieved,  the  surface  of  this  idealized 
ocean  assumes  a  level  surface  of  the  gravity  field.  It  may  be  regarded  as  also 
extending  under  the  continents.  This  level  surface  is  called  the  geoid.  If  the 
potential,  W,  is  given  as  a  function  of  the  coordinates  r,  <{>  and  X,  then  the  geoid  is 
given  by  [Moritz,  1984] 

W  =  W(r,fcX)  =  W0 

The  geoid  is  a  closed  and  continuous  level  surface  which  extends  partially 
inside  the  solid  body  of  the  earth.  The  direction  of  the  gravity  vector  at  any 
point  (plumb  line  or  vertical)  is  normal  to  the  geoid.  The  curvature  of  the  geoid 
displays  discontinuities  at  abrupt  density  variations.  Consequently,  the  geoid  is 
not  an  analytic  surface,  and  therefore  not  a  practical  reference  surface  for 
position  determinations.  The  geoid  however,  is  well  suited  as  a  reference  surface 
for  potential  or  height  differences,  which  are  obtained  by  precise  levelling  in 
combination  with  gravity  measurements. 

To  establish  the  geoid  as  a  reference  surface  for  heights,  the  ocean  water 
level  is  recorded  and  averaged  over  long  intervals  (>  1  year)  using  tide  gauges. 
The  MSL  thus  obtained  represents  an  approximation  to  the  geoid.  The  National 
Geodetic  Vertical  Datum  of  1929  (NGVD  29)  was  derived  for  land  surveys  from 
a  general  adjustment  of  the  first  order  levelling  net  of  both  United  States  of 
America  and  Canada.  In  the  adjustment  MSL  was  observed  at  twenty-one  tide 
stations  in  the  United  States  and  five  in  Canada.  The  geoid  established  by  this 
method  may  deviate  by  ±  1  to  ±  2  m  from  a  level  surface  due  to  periodic, 
nonperiodic  and  secular  variations  [Torge,  1980]. 
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B.  THE  WGS  84  ELLIPSOID 

The  development  of  WGS  84  [DMA,  1987]  was  initiated  by  the  United  States 
Department  of  Defense  for  navigation  and  weapon  systems.  The  Defense 
Mapping  Agency  (DMA)  developed  WGS  84  as  a  replacement  for  WGS  72.  The 
defining  parameters  and  reference  frame  orientation  of  the  WGS  84  Ellipsoid 
and  the  WGS  84  Ellipsoid  Gravity  Formula  are  those  of  the  internationally 
sanctioned  Geodetic  Reference  System  of  1980  (GRS  80).  Accordingly,  a 
geocentric  equipotential  ellipsoid  is  defined  by  the  semimajor  axis  (a),  the 
flattening  (f),  the  equatorial  gravity  (ya),  and  the  angular  velocity  (to).  The  WGS 

84  Ellipsoid  used  the  values: 

a  =  6378137  m 

f  =  1/  298.257223563 

Ya  =  987.03267714  gals 

to  =  7.2921 15xl0'5  rad/s 

The  reference  system  for  GPS  is  WGS  84.  The  precise  geocentric 
coordinates  obtained  from  GPS  receivers  are  in  WGS  84. 

C.  THE  ORTHOMETRIC  HEIGHT,  H,  THE  ELLIPSOID 

HEIGHT,  h  AND  THE  GEOID  HEIGHT,  N 

In  geodetic  applications  three  different  surfaces  or  earth  figures  are 
normally  involved.  First  is  the  earth's  actual  topography;  second,  the  geometric 
surface,  or  ellipsoid  and;  third,  the  equipotential  surface,  the  geoid.  The 
relationship  between  the  earth's  topography,  the  ellipsoid  and  the  geoid  in  a 
section  through  the  earth's  surface  is  shown  in  Figure  1 . 
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Figure  1.  Relationship  between  the  earth's  surface,  the  geoid  and 
the  ellipsoid. 

Two  features  in  the  figure  are  of  particular  interest : 

1.  The  deflection  of  the  vertical,  0,  defined  by  Pizzetti  as  the  angle  at  the 

geoid  between  the  direction  of  the  plumb  line  (normal  to  geoid)  and 
the  normal  to  the  ellipsoid  through  the  point,  P,  on  the  geoid  [Torge, 
1980]. 

2.  The  vertical  separation,  N,  between  the  geoid  and  the  ellipsoid. 

The  deflections  of  the  vertical  and  the  geoid  heights,  N,  depend  on  the 
ellipsoidal  coordinates  and,  hence,  on  the  parameters  of  the  reference  ellipsoid 
and  its  position  with  respect  to  the  earth.  If  they  are  referred  to  the 
geocentrically  situated  mean  earth  ellipsoid,  then  they  are  referred  to  as  absolute 
quantities;  otherwise,  they  are  relative  quantities.  The  absolute  deflections  of  the 
vertical  in  flat  terrain  and  the  highlands  assume  values  between  one  and  ten 
seconds  of  arc;  in  mountainous  areas,  they  vary  between  one-half  and  one  minute 
of  arc.  As  a  result  of  density  variations  within  the  earth,  the  geopotential 
surfaces,  including  the  geoid,  have  irregular  shapes.  Absolute  geoid  heights, 
however,  rarely  exceed  100  m  [Torge,  1980]. 

The  orthometric  heights,  H,  shown  in  Figure  2  are  referred  to  an 
equipotential  surface,  the  geoid.  The  orthometric  height  of  a  point  on  the  earth's 
surface  is  the  distance  from  the  reference  surface  to  the  point,  measured  along 
the  plumb  line  normal  to  the  geoid.  The  ellipsoid  height,  h,  of  a  point  is  the 
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distance  from  the  reference  ellipsoid  to  the  point,  measured  along  the  line  which 
is  normal  to  the  ellipsoid.  For  purposes  of  simplicity,  the  orthometric  height,  H, 
and  the  ellipsoid  height,  h,  are  shown  to  be  along  a  common  vertical.  In  most 
cases,  this  would  cause  a  very  small  error  that  is  considered  insignificant 
compared  to  present  uncertainties  of  the  geoid  height  N  estimates.  The  geoid 
height,  N,  is  defined: 

N  =  h  -  H 


Figure  2.  Relationship  between  the  geoid  height,  N,  the  ellipsoid 
height,  h  and  the  orthometric  height,  H. 


The  orthometric  height  has  greater  physical  meaning  than  the  geometrical 
ellipsoid  height.  The  orthometric  height  has  traditionally  been  determined  by  the 
technique  of  levelling  in  which  increments  of  height  are  obtained  from  the 
intersection  of  the  line  of  sight  of  a  level  instrument,  tangential  to  the 
geopotential  surface  and  passing  through  the  level  axis,  with  two  graduated  rods. 
Accurate  orthometric  height  information  is  needed  for  precise  engineering 
operations  such  as  the  construction  of  dams,  pipelines,  and  tunnels. 

Geoid  heights  have  been  accurately  determined  for  some  major  geodetic 
datums  such  as  the  North  American,  European  and  Australian.  These  datums  are 
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well  supplied  with  astrogeodedc  deflections  and  have  fair  gravity  coverage.  The 
standard  error  of  relative  geoid  heights  in  these  areas  is  about  2  or  3  m. 

D.  THE  GEOID  HEIGHT  FROM  GRAVITY  MEASUREMENT 

Both  the  U.S.  National  Geodetic  Survey  (NGS)  of  National  Oceanic  and 
Atmospheric  Administration  Service  (NOAA)  and  Trimble  Navigation  versions 
of  Rapp’s  360  degree  x  360  order  model  (OSU86F)  were  available  to  me.  Rapp 
computed  two  potential  coefficient  fields  that  are  complete  to  degree  and  order 
360  [Rapp,  1986].  One  field  (OSU86E)  excludes  geophysically  predicted 
anomalies,  while  Rapp's  other  model  (OSU86F)  includes  such  anomalies.  These 
fields  were  computed  using  a  set  of  30-minute  mean  gravity  anomalies  derived 
from  satellite  altimetry  in  the  ocean  areas  and  on  land  from  standard 
measurements. 

Gravity  anomalies  can  be  observed  and  then  used  to  compute  the  geometric 
deviation  of  the  geoid  from  the  ellipsoid.  The  expression  for  the  gravitational 
potential  is  written  in  the  following  form  [Rapp,  1986]: 

kM  r  00  Tall  1  /  \ 

V(r,<j),X)  =  — - -  1  +  X  r  X  (  C  cosmX  +  S  sinmA,  )  P  (sin<J>) 

Y  L  1=2  Lr  J  m=o  V  lm  J  l"1 

where 

r,  (|),  X  :  geocentric  coordinate 

kM  :  geocentric  gravitational  constant 
a  :  equatorial  radius  of  the  reference  ellipsoid 

S^:  fully  normalized  potential  coefficients 

P^  :  fully  normalized  Legendre  function  of  degree  1  and  m 
y  :  normal  gravity 

The  potential  at  a  point,  U,  is  the  scalar  sum  of  V  and  centrifugal  force  potential 
V’  [Ewing,  1976]: 

U(r,«>,X)  =  V(r,4>,X)  + 

The  difference  between  observed  gravity  potential  W  and  the  computed  normal 
gravity  potential  U  is  denoted  by  T  [Moritz,  1984],  so  that 
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W(r,<(>,X)  =  U(r,<j>,X)  +  T(r,<frA) 
compared  the  geoid  defined  by  the  potential  WQ  is  given  by 

W(r,<J>,X)  =  W0 

A  reference  ellipsoid  with  the  same  potential,  WQ  =  UQ,  is  given  by 
U(r,<j),X)  =  WQ 

A  point  P  on  the  geoid  is  projected  onto  the  point  Q  on  the  ellipsoid  along  the 
normal  to  the  ellipsoid  PQ  =  N.  PQ  is  the  distance  between  geoid  and  ellipsoid  at 
the  point  and  is  called  the  geoid  height  (Figure  3). 


U  =  W0 


Figure  3.  Geoid  and  ellipsoid  (Moritz  [1984,  Fig.  2-12]). 
The  disturbing  potential  T  can  be  written  as  [Rapp,  1986]: 


kM  °° 

T(r,<U)  =  -7  1 

7  1=2 


1  i  i  c£,Y,am«iA) 

L1  J  m=0  0=0 
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where 


Q“m=  {Cwa  =  0>  and  5^,0=  1} 

=  {P^sin^osmX,  a  =  0,  and  P^sin^sinmA.,  a  =  1 } 

Since  the  N  is  relatively  small  compared  to  the  ellipsoid  geocentric  radius  R\ 
then  [Ewing,  1976] 

T  =  yN 
and 

T 

N  =  7 

Y 

where  y  is  the  normal  gravity  at  (<j),  A.). 

The  gravity  anomaly,  Ag,  in  the  Molodensky  surface  free-air  anomaly  sense 
is  given  by  [Moritz,  1984]: 

Agp  —  8p  Yq 

The  boundary  condition  that  relates  Ag  and  T  is  [Moritz,  1984] 

_  3T  I_  8Tt 

®  8h  r  8h 

where  h  is  the  distance  along  the  plumb  line  direction.  Neglecting  deflections  of 
the  vertical  we  have  [Rapp,  1986] 


^  =-7t-£o- n^T  i  i < cs,- - ea > y^*.x) 

Y  1=2  Lr-I  m=0  o=0 


where  ^  ( i  =  h,  y)  are  ellipsoidal  corrections. 

The  Stokes  function,  S(\j/),  can  then  be  used  to  solve  the  geoid  heights  above  the 
geodetic  ellipsoid  [Ewing,  1976], 
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Ng=  4^V'/o7olAg(V,a)S(v)sin¥dvdCt 

*  m 

where 

Ym  :  the  mean  value  of  normal  gravity 

y  :  the  angular  distance  between  the  point  where  N  is  being 
determined  and  the  area  where  the  effect  of  Ag  is  being 
considered 

a  :  the  azimuth  from  the  affected  point  to  that  causing  the  effect 
Ag  :  the  gravity  anomaly 
S(y)  :  the  Stokes  function 

The  Stokes  function  is  given  by  [Ewing,  1976] 

S(y)  =  esc  2^  +  1  -  6  sin  ^  -  5  cosy  -  3  cosy  ln(sin  +  sin2 

The  gravitational  potential  using  degree  360  and  order  360  has  been  tested 
through  comparison  of  Doppler  station  geoid  heights  with  geoid  heights  from 
Rapp's  versions  of  geopotential  models.  The  agreement  between  the  two  geoid 
height  measurements  is  approximately  ±  1.6  m  [Rapp,  1986]. 
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A.  THE  GLOBAL  POSITIONING  SYSTEM 

The  Global  Positioning  System  (GPS)  calls  for  a  precise  navigation  system 
divided  into  three  segments:  space  segment,  control  segment  and  user  equipment. 
The  space  segment  will  consist  of  three  orbital  planes  of  satellites  at  inclinations 
of  120°  in  circular  orbits  at  altitudes  of  20,000  km.  Each  plane  will  eventually 
contain  six  to  eight  satellites  to  give  the  three  dimensions  of  position,  velocity  and 
precise  time  24  hours  a  day  anywhere  in  the  world.  The  control  segment  consists 
of  the  ground  stations  necessary  to  track  the  satellites,  monitor  the  system 
operation,  and  periodically  provide  corrections  to  the  navigation  and  time 
signals.  Each  satellite  broadcasts  signals  containing  information  on  its  position. 
The  GPS  satellite  nansmits  signals  at  two  L-band  frequencies  (1227  and  1575 
MHZ)  to  permit  corrections  for  ionospheric  corrections.  The  signals  are 
modulated  with  two  codes:  P,  which  provides  for  precise  measurement,  and  C/A, 
which  permits  easy  lock-on  to  the  desired  signal.  The  user  segment  consists  of 
the  equipment  necessary  to  convert  the  satellite  navigation  message  into  useful 
navigation  information. 

The  navigation  message  contains  the  data  that  the  user’s  receiver  requires  to 
perform  the  operations  and  computations  for  successful  navigation  with  GPS. 
The  data  includes:  (1)  information  on  the  status  of  the  Space  Vehicle  (SV);  (2) 
time  synchronization  information  for  the  conversion  of  the  C/A  to  P  code;  (3) 
parameters  for  computing  the  clock  correction  and  the  ephemeris  of  the  SV;  and 
(4)  corrections  for  delays  in  the  propagation  of  the  signal  through  the 
atmosphere.  In  addition  the  data  contain  almanac  information  to  define  the 
approximate  ephemeris  and  to  give  the  status  of  all  the  other  SV  information 
which  is  required  for  use  in  signal  acquisition  [Milliken,  1980].  Ranges  to  the 
satellites  are  determined  by  signal  transit  times  multiplied  by  the  speed  of  light 
(299,792,458  m/sec).  The  transmitted  message  contains  ephemeris  parameters 
that  enable  the  user  to  calculate  the  position  of  each  satellite  at  the  time  of  the 
transmission  of  the  signal. 


B.  CARRIER  PHASE  MEASUREMENTS 

GPS  measurements  can  be  made  using  the  pseudo-range  and  the  carrier 
phase.  The  pseudo-range  is  essentially  a  measurement  of  distance  contaminated 
by  clock  error.  When  four  satellites  are  observed  simultaneously,  the  three 
dimensional  position  of  the  ground  receiver  can  be  determined  along  with  the 
receiver  clock  offset  at  a  single  epoch.  The  accuracy  of  pseudo-ranges  is  affected 
by  multipath  effects,  which  depend  on  the  antenna  design,  and  its  height  above  the 
ground. 

Carrier  phase  measurements  are  more  precise  than  pseudo-range  and  are 
not  as  vulnerable  to  multipath  effects.  They  can  be  used  to  compute  the  precise 
base  line  components  AX,  AY,  AZ  between  two  receivers.  Phase  measurements 
are  made  by  beating  the  received  carrier  with  the  signal  from  a  local  oscillator 
internal  to  the  GPS  receiver.  The  slant  range  from  a  GPS  receiver  to  a  satellite 
can  be  modelled  in  terms  of  time.  It  takes  the  signal  time  to  travel  between  the 
satellite  and  the  receiver  or  the  equivalent  number  of  cycles.  The  range  of  the 
cycles  will  consist  of  an  integer  and  fractional  number  of  cycles.  When  a 
receiver  locks  onto  the  carrier  signal,  it  can  immediately  measure  the  fractional 
part  and  begin  counting  subsequent  integer  cycles,  but  it  can  not  measure  or 
account  for  the  initial  integer  number  of  cycles  that  preceded  the  initial  fractional 
part.  The  initial  integer  ambiguity  which  biases  the  subsequent  measurements  is 
called  the  initial  integer  ambiguity  bias. 

GPS  uses  a  one-way  carrier  beat  phase.  The  GPS  satellite  and  receiver  are 
controlled  by  separate  clocks.  The  satellite  clock  generates  the  signal,  and  the 
receiver  clock  detects  when  the  signal  arrives.  An  error  in  the  synchronization 
of  the  clocks  of  1  microsecond  creates  an  error  in  range  of  300  m. 

C.  ONE-WAY  CARRIER  PHASE  MEASUREMENT 
DIFFERENCING 

A  single  difference,  SD(j,i),  is  formed  by  differencing  carrier  beat  phase 
observables  from  two  receivers  1  and  2  at  the  same  observation  epochs  i  of  same 
satellite  j.  The  equation  is  given  by  [Remondi,  1984]: 

SD(j,i)  =  S(2,j,i)  -  S(l,j,i) 
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where  S(kj,i)  is  the  raw,  unpreprocessed,  fractional  phase  plus  the  count  made  at 
epoch  i  by  receiver  k  for  satellite  j.  The  main  advantage  of  the  single  difference 
is  that  it  reduces  or  eliminates  satellite  orbital  and  clock  errors,  because  they  are 
common  to  both  receivers.  Its  disadvantage  is  that  one  can  not  exploit  the  integer 
nature  of  integer  ambiguities.  Thus,  for  short  base  lines,  the  ultimate  in  accuracy 
may  not  be  achievable  [Remondi,1985]. 

A  double  difference,  DD(j,k,i),  is  formed  by  differencing  single  differences 
between  a  reference  satellite  j  and  another  satellite  k  at  the  same  epoch  i.  The 
equation  is  given  by  [Remondi,  1984]: 

DD(j,k,i)  =  SD(k,i)  -  SD(j,i) 

The  advantage  of  the  double  differences  is  that  the  receiver  clock  dependent 
terms  are  eliminated  because  the  differences  for  each  epoch  are  correlated.  The 
significance  of  the  removal  of  those  terms  is  to  reduce  from  nanoseconds  to 
microseconds  the  timing  accuracy  required  to  achieve  one  cycle  accuracy.  For 
short  base  lines,  the  integer  ambiguities  can  be  isolated,  since  the  contribution 
made  by  the  clock  drift  is  reduced  [Remondi,  1985].  The  Trimble  4000SX 
receiver  achieves  sub-microsecond  accuracy  by  using  the  C/A  code  timing 
information  [Ashjaee,  1985]. 

A  triple  difference,  TD(j,k,i),  is  formed  by  differencing  the  double 
differences  for  the  same  satellite  pair  at  some  integer  of  succeeding  epochs  i+1 
[Remondi,  1984]. 

TD(j,k,i)  =  DD(j,k,i+l)  -  DD(j,k,i) 

The  advantage  of  the  triple  difference  is  that  it  eliminates  all  the  time  independent 
terms,  namely  the  initial  integer  ambiguities,  and  becomes  insensitive  to  the 
initial  ambiguities  and  any  cycle  slips  when  the  receiver  loses  lock.  The 
disadvantage  of  the  triple  difference  is,  another  level  of  correlation,  loss  of 
resolution  and  a  greater  number  of  observations.  Triple  differences  are  already 
correlated  with  respect  to  satellite  because  of  the  underlying  double  differences 
and  are  further  correlated  with  respect  to  time  because  consecutive  triple 
observations  will  have  common  DD(j,k,i+l)  terms. 
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For  short  base  lines  local  area  integer  ambiguities  can  easily  be  resolved 
because  unmodelled  errors  are  highly  correlated  between  the  two  antenna  sites  and 
are  mostly  eliminated  by  differencing.  Algorithms  can  take  advantage  of  the 
integer  nature  of  the  initial  ambiguities  and  solve  for  them  [Remondi,  1984]. 

D.  GEOID  HEIGHT  FROM  GPS 

Surface  Fitting  techniques  can  be  used  with  GPS-derived  geoid  heights.  GPS 
stations  are  likely  to  be  close  together,  of  the  order  of  a  few  tens  of  km,  and  the  local 
geoid  can  be  estimated  directly  if  levelling  data  is  available  in  the  area.  This 
together  with  it’s  relative  simplicity  makes  the  method  practical  for  correcting  GPS 
heights. 

It  is  possible  to  use  the  two  sets  of  elevations  (that  is,  levelled  and  GPS- 
derived)  to  define  two  distinct  planes.  The  published  levelled  elevations  are 
referred  to  the  geoid  and  the  GPS  elevations  are  referred  to  an  ellipsoidal  surface. 
If  both  elevations  are  made  equal  at  one  bench  mark  (control  station),  then  in 
general,  the  other  bench  marks  will  have  two  elevation  values.  After  several 
different  models  were  studied  two  mathematical  surface  models  were  chosen  in  this 
study.  The  five-parameter  model  found  suitable  for  use  on  the  NPS  campus  is 

N0(X,Y,Z)  =  h0  +  Ah  -  H  +  aAY  +  bAX2  +  cAY2  +  dAXAY 

and  the  six-parameter  model  selected  for  a  larger  area  near  Monterey,  California,  is 
given  by 


R 


N0(X,Y,Z)  =  h0  +  Ah  -  H  +  a’ AX  +  b'AZ  +  c’AX2  +  d’AY2  +  e’AXAY 


where 

N0(X,Y,Z) :  the  global  geoid  height,  obtained  by  gravity  measurements 
h0  :  the  ellipsoid  height  of  the  reference  point,  including  a  constant 
correction  to  N0(X,Y ,Z)  at  the  reference  station 
Ah  :  the  ellipsoid  height  difference  with  respect  to  the  reference  station 
H  :  the  published  or  levelled  elevations  at  each  control  station 
a,  b,  c,  d,  a’,  b\  c’,  d\  e’ :  the  coefficients  to  be  determined 
AX,  AY,  AZ  :  the  coordinate  differences  in  WGS  84 

R 
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These  models  can  be  solved  to  determine  the  east-west  and  north-south  tilts 
which  are  absorbed  by  two  rotations,  one  around  the  north  axis  (AY)  and  the 
other  around  the  east  axis  (AX)  in  the  horizon  system.  Their  separation  is 
absorbed  by  the  scale  correction  [Zilkoski,  1988].  The  local  geoid  heights,  N, 
can  be  found  from  the  equation 

N  =  N0(X,Y,Z)  +  AN 

where  AN  is  the  variation  of  geoid  height  in  local  area.  The  equation  is  given  by 

AN  =  H  +  N0(X,Y,Z)  -  Ah 

AN  =  h0+  aAY  +  bAX2  +  cAY2  +  dAXAY 
or 

AN  =  h0  +  a’AX  +  b'AZ  +  c’AX2  +  d'AY2  +  e'AXAY 

To  solve  for  h0,  a,  b,  c,  d,  a',  b',  c\  d’  and  e’,  the  global  geoid  heights 
N0(X,Y,Z)  can  be  obtained  from  the  Geoid.exe  program  described  in  Chapter 
IV.  The  geoid  model  can  be  rearranged  to  give 

H  =  h0  +  Ah  -  N0(X,Y,Z)  +  aAY  +  bAX2  +  cAY2  +  dAXAY 
or 

H  =  h0  +  Ah  -  N0(X,Y,Z)  +  a'AX  +  b'AZ  +  c’AX2  +  d'AY2  +  e’AXAY 

Then  h0  a,  b,  c,  d,  a',  b\  c\  d'  and  e'  can  be  solved  by  the  least  squares  method. 
The  geoid  height,  N,  found  by  this  method  appears  to  be  adequate  for  areas  up  to 
50  km  x  50  km  where  the  geoid  is  smooth  [King,  1985]. 

E.  METHOD  OF  COLLOCATION 

The  method  of  collocation  was  derived  from  least  squares  interpolation. 
This  method  was  used  to  predict  the  geoid  height,  N,  for  a  local  area.  The  geoid 
model  is  given  by 

N  =  N0(X,Y,Z)  +  AN  +  S  +  n 
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where 

S  :  the  signal 
n  :  the  noise 

N0(X,Y,Z)  +  AN  :  the  system  function  from  the  surface  models  at  control 
points  where  both  h  and  H  are  known 

Then,  for  a  given  N0(X,Y,Z)  the  method  of  collocation  can  be  used  to  determine  Sq 
at  these  control  points  and  to  predict  Sp  at  other  points  in  the  local  area.  At  any  point 
in  the  local  area  the  value  of  h  can  be  determined  by  using  differential  GPS 
measurements  between  a  known  point  and  any  other  point  by  the  equation 

h  =  h0  +  Ah 

The  value  of  H  at  any  point  is  given  by 

H  =  h  -  N0(X,Y,Z)  -  AN  -  S  -  n 

The  value  of  H  at  any  point  in  the  local  area,  which  depends  on  the  value  of  h  at 
the  control  points  and  the  differential  GPS  measurements,  can  be  predicted  with  an 
accuracy  of  ±  n.  The  general  form  of  the  observation  equation  in  the  method  of 
collocation  is  [Jeyapalan,  1977]: 

x  =  A«X  +  S  +  n  +  0*S 
q  q  p 

where 

x  :  the  vector  of  the  observation  (x  =  Ah  -  N0(X,Y,Z)  -  H) 

A  :  a  given  rectangular  coefficient 

X  :  the  vector  of  the  systematic  parameters  (h0,  a,  b,  w ,  d,  or  h0,  a',  b',  c', 
d\  e') 

Sq  :  a  signal  vector  at  q  observation  points 
nq  :  a  vector  of  measuring  errors,  noise  at  q  points 
Sp  :  a  signal  vector  at  p  unknown  stations 
•  :  indicates  matrix  multiplication 
0  :  the  null  matrix 
If 


•i 


■ 

i 


« 


l 


Z  —  S  +  n 
q  q  q 
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Then 

x  —  A*X  +  Z  +  0*S 
q  q  p 

then 

X  =  (ATCq1A)‘1  ATCq'1x 

and 

S  =  C  C  1  (x  -  A  X) 
p  pq  q 

where  the  variance-covariance  matrix  C  of  the  Z  and  S  vectors  is 

q  p 


The  essence  of  this  method  is  that  by  some  means  a  covariance  matrix  can 
be  assigned  to  the  signal.  For  noise  it  will  be  possible  to  assign  a  diagonal  weight 
matrix. 

Sp  are  the  values  of  the  signal  at  the  interpolated  stations.  Suppose  there  are 
q  observations  (and  values  of  Sq),  p  interpolated  values  of  Sp  and  m  model 
parameters.  The  covariance  matrices  are  the  following  [Bomford,  1980]: 

(i)  Cq,  the  expected  covariance  between  the  observed  x’s  for  all  pairs  of 

the  q  observations.  It  is  a  q  x  q  matrix. 

(ii)  C„,  the  expected  covariance  between  the  signals  for  all  pairs  of  the  q 

observations.  It  is  a  q  x  q  matrix. 

(iii)  Cnq,  the  expected  value  of  n  at  each  station.  It  is  a  q  x  q  diagonal 

matrix. 

(iv)  C  ,  the  expected  covariance  between  all  pairs  of  mixed  observed 

and  interpolated  signals.  It  is  a  p  x  q  matrix. 

(v)  C  ,  the  expected  covariance  between  the  signals  at  pairs  of 

interpolated  stations.  It  is  a  p  x  p  matrix. 

The  variance  of  the  noise,  Cnqcan  be  estimated  in  the  usual  way  according  to 
the  circumstances  at  each  station,  different  types  of  instrument,  etc.  The  noise  at 
different  points  is  independent;  hence 
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where  on  is  the  standard  error  of  the  noise. 

The  expected  covariance  between  the  observation  x’s,  Cq  can  be  obtained 
by 


because  the  signal  is  small. 

The  can  be  be  computed  from  a  simple  function  whose  parameters  can 
be  determined  by  using  empirical  data.  The  covariance  function  is  positive  and 
definite.  These  two  characteristics  are  found  in  many  functions.  Functions 
commonly  used  for  covariance  may  be  constant,  sinusoidal,  Gaussian, 
exponential,  exponential  cosine,  and  exponential  sine  and  cosine.  In  this  thesis  the 
sinusoidal  function  was  used.  The  function  is  given  by 

C  =  B  sin(kr) 

sq 

C(r)  =  C  +  B  sin(kr) 

where  Cnq  is  the  standard  deviation  of  the  control  stations  and  B,  k  are 
coefficients  to  be  determined. 

The  parameters  can  be  determined  by  least  squares,  and  residuals  at  each 
point  can  be  computed.  From  the  residual,  the  covariance  between  points  at  a 
distance  (r)  can  be  computed  by 


C(r)  = 


y  v  v 

La  0  r 


i=l _ 

n  -  1 


where  V0  is  the  residual  at  the  center,  Vr  is  the  residual  at  a  point  which  is  at  a 
distance  r  from  the  center  and  n  is  the  number  of  points.  There  are  several 
methods  of  determining  C(r)  of  which  the  concentric  circle  approach  is  the 
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simplest.  The  coefficients  B  and  k  can  then  be  computed  from  the  computed 
covariance. 

The  expected  covariance  between  the  signals  at  pairs  of  interpolated  stations 
Cp  can  be  obtained  by  substituting  the  distance  between  the  interpolated  stations 

and  control  stations  into  the  covariance  function. 

The  expected  covariance  between  all  points  of  mixed  observed  and 
interpolated  signals  CM  can  be  obtained  from  the  covariance  function  for  each 

distance  from  the  interpolated  stations  to  the  control  stations. 

In  this  thesis  two  cases  were  studied  in  solving  for  the  geoid  height  model. 

Case  1.  Cq  =  I  and  =  0 

Case  2.  Cq  *  I  and  C^  *  0 

I  start  with  assuming 

C  =  P 1  =  I 
q 


where  P  is  the  weight  matrix  and  I  is  the  unit  matrix. 

Then 

X0=(AtP  A)1  ATP  x 
and 

S  =  0«P  (x  -  A  Xn)  =  0 

p  0 

Using  X  to  compute  residuals,  v,  and  then  estimating  C  .  C  and  P,  it  is  then 

q  pq 

found  that 

X  =(AtC1A)1  AtC  1  x 

1  q  q 

and 

S  =  C  C  1  (x-  AX.) 
p  pq  q  i 
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IV. 


A.  PLANNING 

Nine  temporary  bench  marks  were  established  on  the  NPS  campus  (Figure  4). 
Bench  marks  GH7  and  GH8,  designated  as  check  marks,  were  established  in  the 
center  of  the  levelling  loop.  H  of  bench  marks  was  obtained  by  differential 
levelling.  Ah  was  measured  with  a  pair  of  GPS  receivers. 
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Eight  permanent  bench  marks  were  recovered  in  the  study  area  (Figure  5) 
(Vertical  Control  Data,  1961].  Bench  mark  GWM  27,  designated  as  a  check  mark, 
was  roughly  in  the  center  of  the  survey  area  for  the  permanent  bench  marks.  H  was 
obtained  from  the  published  elevations.  Ah  was  also  measured  with  a  pair  of  GPS 
receivers. 
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B.  DATA  COLLECTION 
1.  Obtaining  H 

a.  Precise  Levelling 

Precise  levelling  was  run  on  the  NPS  campus  on  February  19,  1988, 
and  March  24,  1988.  The  elevation  of  station  TREE  was  assumed  to  be  zero 
above  the  MSL,  so  the  orthometric  height  differences,  AH,  for  each  temporary 
bench  mark  could  be  observed.  The  differences  of  the  relative  elevations,  AH, 
for  the  forward  and  the  backward  sights  on  the  NPS  campus  are  listed  in  Table  1 . 


Table  1.  AH  FOR  THE  NPS  CAMPUS 


Bench 

mark 

AH  (m) 

From 

To 

Forward 

Backward 

Mean 

TREE 

GH1 

-0.808 

0.808 

-0.808 

GH1 

GH2 

5.388 

-5.384 

5.386 

GH2 

GH3 

3.150 

-3.149 

3.149 

GH3 

GH4 

1.535 

-1.536 

1.535 

GH4 

GH5 

-1.016 

1.016 

-1.016 

GH5 

GH6 

-4.214 

4.215 

-4.214 

GH6 

GH7 

0.489 

-0.489 

0.489 

GH7 

TREE 

-4.031 

4.033 

-4.032 

GH6 

GH2 

0.547 

-0.546 

0.547 

GH7 

GH2 

0.057 

-0.057 

0.057 

GIB 

GH8 

0.716 

-0.718 

0.717 

GH8 

GH4 

0.817 

0.817 

0.817 

To  avoid  an  obstruction  near  bench  mark  D  697,  which  is  10  ft  north 
of  a  12-inch  cedar  tree,  the  temporary  bench  mark  TEMP  1  was  established 
nearby.  Offset  levelling  was  run  on  April  8,  1988.  The  elevation  of  TEMP  1 
was  offset  by  -0.401  m  from  the  elevation  of  D  697. 

F  813  which  was  set  in  the  west  end  of  the  south  abutment  of  a  steel 
bridge  over  the  Salinas  River  was  offset  to  temporary  bench  mark  TEMP  2.  The 
offset  levelling  was  run  on  April  9,  1988.  The  elevation  of  TEMP  2  was  offset 
by  -2.218  m  from  the  elevation  of  F  813. 
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R 

b.  Published  Elevation 

H  values  for  the  larger  of  our  two  study  areas  in  and  near  Monterey, 
i  California,  areas  were  obtained  from  the  published  elevations  printed  by  U.S. 

Department  of  Commerce  Coast  and  Geodetic  Survey  Washington  D.C.  [1961]. 
The  geodetic  datum  used  in  this  publication  was  the  NGVD  29.  First-order  spirit 
levelling  has  extended  this  datum  over  most  of  the  continent.  Although  first- 
I  order  lines  may  be  300  km  apart  in  some  western  areas,  most  points  in  the 

country  are  no  more  than  50  km  from  an  estimated  first-order  bench  mark.  A 
readjustment  of  the  this  network  is  the  NAVD  88.  This  new  adjustment  will  be 
made  to  the  geopotential  surface  rather  than  to  sea  level,  and  it  will  place  the 
I  existing  vertical  data  in  a  form  that  makes  it  most  consistent  and  accessible  to  the 

user.  Changes  to  older  published  elevations  are  not  expected  to  exceed  15 
decimeters  [NASA,  1978].  Table  2  lists  the  orthometric  heights  of  the  permanent 
bench  marks  we  used. 

f  Table  2.  H  FOR  PERMANENT  BENCH  MARKS  UESD  IN  THIS 

|  STUDY 

li 

\ 
i 

i 

i 


Bench  mark 

Hftn) 

K  152 

16.965 

B  21 

5.787 

S  812 

15.978 

P  812 

13.134 

D  697 

30.057 

J  697 

85.061 

F  813 

7.983 

GWM27 

63.487 

2.  Obtaining  Ah 

a.  Satellite  Observation  Plan 

Due  to  the  positions  of  the  satellite  orbits  during  this  study,  the 
observing  window  was  between  2100  and  0200  hours  Pacific  Standard  Time 
(PST=UTC+8  hours.).  The  time  period  was  between  February  2,  1988,  and  May 
12,  1988. 
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b.  Satellite  selection 

The  same  five  satellites  (SV)  (6,  9,  11,  12  and  13)  were  used  for  all 
observations. 

c.  Position  Dilution  of  Precision  (PDOP) 

The  accuracy  to  which  positions  are  determined  using  GPS  depends 
on  two  factors:  (i)  satellite  configuration  geometry,  and  (ii)  measurement 
accuracy.  GPS  measurement  accuracy  represents  the  combined  effect  of 
ephemeris  uncertainties,  propagation  errors,  clock  and  timing  errors,  and 
receiver  noise. 

The  effect  of  satellite  configuration  geometry  is  expressed  by  the 
dilution  of  precision  (DOP)  factor,  which  is  the  ratio  of  the  positioning  accuracy 
to  the  measurement  accuracy  [Wells,  1987]. 

a  =  DOP  •  o0 


where 

ct0  :  is  the  measurement  accuracy  (standard  deviation) 
a  :  is  the  positioning  accuracy  (standard  deviation  in  one  coordinate) 
The  value  of  GDOP  itself  is  a  composite  measure  that  reflects  the 
influence  of  satellite  geometry  on  the  combined  accuracy  of  the  estimation  of 
observation  time  (user  clock  offset)  and  receiver  position  [Milliken,  1986]. 

GDOP2=  PDOP2  +  TDOP2 

TDOP  is  the  Time  Dilution  Of  Precision,  the  error  in  the  clock  of  the  receiver 
bias  multiplied  by  the  velocity  of  light.  The  four  best  satellites  selected  by  the 
receiver  are  those  with  the  lowest  GDOP.  Trimble  recommends  that  the  rapidly 
changing  PDOP  provides  better  geometry  for  phase  differencing  techniques. 
Low  or  constant  PDOP  provides  weaker  solutions.  The  4000SX  receiver  does 
not  record  GDOP,  but  it  does  record  PDOP  every  five  minutes.  PDOP  was 
about  4.5  for  all  GPS  observations. 
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C.  INSTRUMENTATION 

1.  Levelling 

A  Zeiss  model  Ni-2  level  (Ser.  #  82377)  was  used  at  the  temporary 
bench  marks  for  the  third-order,  class  I  levelling.  It  has  a  32-diameter 
magnification,  produces  an  erect  image,  and  has  stadia  constants  of  333  or  100. 
A  Peg  test  was  performed  before  levelling.  The  level  error,  or  c-value  was  also 
checked  before  the  beginning  the  levelling.  The  c-value  was  -0.006  mm/m  which 
was  less  than  +0.05  mm/m,  so  it  was  not  necessary  to  adjust  the  level  [Bodnar, 
1975].  The  level  contains  a  bubble  tube  to  permit  positioning  parallel  to  the 
geoid.  When  properly  set  up  at  a  point,  the  telescope  is  locked  so  that  it  will 
rotate  through  a  360°  arc  in  a  horizontal  plane.  With  the  level  locked  in  position 
readings  are  made  on  two  calibrated  staffs  held  in  upright  positions  ahead  of  and 
behind  the  instrument.  The  difference  between  readings  is  the  difference  in 
elevation  between  points.  Dietzgen  Model  #  6450  metric  rods  were  used  in  the 
levelling.  These  rods  are  graduated  in  centimeters.  The  actual  reading  is 
estimated  to  the  nearest  1  mm.  Rod  levels  were  used  to  indicate  when  were 
vertical. 

2.  GPS 

a.  Trimble  4000SX  Receiver 

A  complete  description  of  the  Trimble  4000SX  receiver  is  given  by 
Trimble  Navigation  [Trimble,  1987a].  NPS  operates  three  Trimble  4000SX  GPS 
receivers.  For  this  study  I  fixed  one  antenna  to  the  roof  of  Building  224  on  the 
NPS  campus,  and  one  was  carried  to  the  field.  The  4000SX  is  capable  of 
observing  the  C/A  code,  integrated  Doppler,  and  carrier  beat  phases  of  up  to  five 
satellites  simultaneously.  Its  ability  to  use  the  C/A  code  allows  the  receiver  to  be 
used  as  a  stand  alone  navigation  system  to  determine  positions  using  Doppler- 
smoothed  pseudoranges  and  velocities  [Ashjaee,  1985],  The  receiver  uses  the 
C/A  code  in  a  time  transfer  mode  to  determine  the  offset  and  drift  of  its  own 
clock  and  thus  provide  accurate  time  tags  for  the  observations  without  the 
requirement  of  an  external  atomic  clock  or  synchronization  with  the  receiver  at 
the  other  end  of  the  baseline.  The  reference  position  (the  geodetic  coordinates  of 
the  antenna)  and  the  practical  options  chosen  must  be  entered  into  the  receiver  via 
the  receiver  key  pad.  The  4000SX  requires  115  V  AC  power  or  12  V  DC 
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power  supply.  115  V  AC  power  was  used  for  the  NPS  campus  measurements. 
12  V  EXT,  supplied  by  a  car  battery,  was  used  in  the  field. 

b.  Grid  Personal  Computer  (PC) 

For  precise  relative  positioning  the  4000SX  receiver  transmits  data 
through  an  RS-232  port  to  a  microcomputer  (Grid  PC)  for  storage  on  floppy 
disks  for  post-processing.  The  Grid  PC  uses  either  115  V  AC  or  12  V  DC.  On 
the  NPS  campus  a  115  V  AC  power  supply  was  used,  so  measurements  could  be 
made  for  four  hours.  The  12-volt  battery  in  the  Grid  PC  lasts  about  120 
minutes,  so  measurements  were  taken  for  only  100  minutes  in  the  field. 

c.  Antennas 

Multipath-resistant  Trimble  microstrip  antennas  were  installed  on 
Building  224  on  the  NPS  campus  and  over  the  various  bench  marks.  The  antenna 
heights  from  the  center  of  bench  marks  to  the  edge  of  the  antenna's  ground  plane 
were  measured  before  and  after  GPS  observations.  The  field  antenna  was 
mounted  on  a  tripod  with  a  tribach  and  optical  plummet  for  centering  and 
levelling  the  antenna.  Arrows  on  the  antennas  were  oriented  to  the  north  at  both 
stations  using  a  magnetic  compass. 

d.  Meteorological  Instruments 

A  barometer  and  a  sling  psychrometer  were  used  to  measure 
atmospheric  pressure,  relative  humidity  and  air  temperature  at  each  field  site. 

e.  The  steel  tape 

A  three-meter  steel  tape  was  used  to  measure  the  antenna  height. 
This  tape  can  be  read  to  1  cm.  Readings  are  estimated  to  1  mm. 

D.  SOFTWARE 

A  complete  description  of  the  Trimble-supplied  Trimvec  software  is  to  be 
found  in  Trimble  Navigation  [Trimble,  1987b].  The  software  provides  data 
logging,  baseline  computation  and  datum  transformation  programs.  These 
operate  with  an  IBM  compatible  personal  computers.  A  printer  and  a  hard  disk 
are  used  with  the  planning  and  processing  programs.  The  following  programs 
are  used  on  3-1/2  inch  micro-floppies: 

1.  Data  Logger  on  Disk  1 

In  relative  positioning,  the  receiver  is  controlled  from  the  Grid  PC  by 
version  D  of  Trimble's  Gridlog5.bat  program.  Each  observation  session  was 
initialized  to  log  data  when  a  minimum  of  four  satellites  were  15°  above  the 


antenna's  horizon.  Five  satellites  were  designated  for  each  observing  session. 
The  observables  and  receiver  clock  parameters  were  logged  to  a  floppy  disk 
every  15  seconds  and  the  C/A  code-determined  antenna  position  and  PDOP  every 
five  minutes.  The  GPS  navigation  message  was  logged  to  a  separate  file  at  the 
beginning  of  the  session. 

2.  Post-processor  on  Disk  2 

The  data  was  processed  using  the  Trimvec  Trim640  program.  Revision 
AB.  Trim640  is  a  relative  positioning,  post- processing  program  that  provides 
triple  and  double  difference  solutions  for  two  sets  of  the  4000SX  carrier  phase 
data  logged  simultaneously  at  two  stations.  Trim640  adopts  the  best  C/A  code 
position  during  the  data  loading.  Only  the  broadcast  ephemeris  can  be  used  to 
compute  fixed  orbit  satellite  positions.  Trim640  limits  processing  to  700 
epoches,  so  the  first  700  epoches  for  each  observing  session  are  used.  The 
antenna  on  the  roof  of  Building  224  on  the  NPS  campus  was  used  as  reference 
station  and  its  coordinates  were  kept  fixed.  Permanent  bench  marks  were  chosen 
as  Trim640  reference  stations  for  the  off  campus  area  because  the  observing 
sessions  at  them  were  shorter  than  the  observing  session  at  Building  224  on  the 
NPS  campus. 

3.  Geoid.exe  on  Disk  3 

This  program  computes  global  geoid  height,  N0(X,Y,Z),  at  any  point 
with  an  accuracy  of  a  few  meters  [Trimble,  1987b].  The  global  geoid  heights, 
N0(X,Y,Z),  are  based  on  Rapp's  1978  360  x  360  model,  an  harmonic  expansion 

of  the  geopotential  referred  to  GRS  80.  Heights  are  suitable  for  showing  relative 
shape  and  trends  in  geoid.  Global  geoid  height  for  the  Monterey  Bay  area  from 
36°  34'  N  to  36°  49'  N  latitude,  and  from  121°  34'  W  to  121°  55’  W  longitude  is 
given  in  Figure  6.  A  Fortran  program  GPSCON  (Appendix  A)  was  used  for  the 
contoured  plot. 

A  second  program  also  named  Geoid.exe  was  obtained  from  the  NGS  of 
NOAA’s  National  Ocean  Service  Charting  and  Geodetic  Services  uses  the  same 
Rapp’s  1986  360  x  360  model.  The  global  geoid  height  in  the  Monterey  Bay  area 
from  the  NGS  program  is  shown  in  Figure  7.  Because  the  NGS  program  rounds 
off  to  the  nearest  decimeters,  there  are  some  steps  on  the  curves.  The  differences 
between  Figure  6  and  7  may  be  attributed  to  differences  between  Rapp’s  1978 
and  1986  models. 
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LONG ITUOE  121  OEGREE  WEST  (MINUTE) 


Figure  6.  The  global  geoid  height  in  the  Monterey  Bay  area  calculated 
using  the  Trimvec  Geoid.exe  program. 


Figure  7.  The  global  geoid  height  in  the  Monterey  Bay  area  calculated 
using  the  NGS  Geoid.exe  program. 


The  program  Stvis.exe  is  a  system  for  generating  visibility  data  for  the 
GPS  satellites.  The  program  takes  as  input  the  observer's  latitude,  longitude,  a 
mnemonic  name  for  the  location,  data  for  calculations  and  optionally  the  time  zone 
offset  from  Universal  Time  Coordinates  (UTC).  NPS  in  Monterey,  California,  was 
used  as  the  reference  station  (  36°  35'  56.1"  N,  121°  52'  36.0"  W,  and  altitude  -19 
m).  Satellite  orbit  data  is  contained  in  the  file,  Almanac.dat,  which  is  produced  by 
processing  data  collected  from  the  actual  satellites.  The  satellite  visibility  chart 
shows  the  tracks  of  the  GPS  satellites  during  the  planned  observation  session 
(Figure  8). 


Figure  8.  Sky  plots  of  satellite  tracks  for  the  Monterey  Bay  area. 
Elevation  angles  are  dotted  concentric  circles.  Zenith  is  at  the  center. 

29 


A.  ORTHOMETRIC  HEIGHT  ON  NPS  CAMPUS 

In  calculating  the  most  probable  elevations  for  the  temporary  bench  marks, 
station  TREE  was  used  as  a  reference,  where  the  elevation  was  assumed  to  be  0 
meters  above  MSL.  The  interlocking  levelling  circuit  on  the  NPS  campus  is 
shown  in  Figure  9. 


Figure  9.  The  interlocking  levelling  circuit  on  the  NPS  campus. 


A  constrained  least  squares  adjustment  was  used  to  do  the  computation.  A 
weight  of  100  was  assumed  at  the  TREE  and  weights  of  1  were  used  for  all 
observations.  The  13  observation  equations  are: 
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P,  Htrhb  -  P,  <  0.000  +  V,  ) 

P2<H,-»W  =  P2  ( -0.808  +  V2  ) 

p3  (  Ha  -  H,  )  =  P3  (  5.386  +  V,  ) 

p4  (  h3  -  H2  )  =  P4  (  3.149  +  V4  ) 

P5  (  H4  -  H3  )  =  P5  (  1.535  +  V5  ) 

V  H5  -  H4  )  =  P6(  -1.016  -  v6  ) 

P,  (  H6  -  Hs  )  =  P?  (  -4.214  +  V,  ) 

P,  <  H,  '  H6  >  =  P*  <  0489  +  V,  > 

P,<tWH7  >  =  P,  (  -4.032  +  V,  , 

P,o<  -  H,  )  -  P,„(  0.547  +  V1Q) 

pn(  Hj  -  H,  )  =  P„(  0.057  +  V„) 

P,/  H8  -  H3  )  =  P,2(  0.717  +  V,2) 

P13<  H4  '  H«  >  =  P,3<  O-817  +  V,3> 

The  observation  equations  in  matrix  form  are  [Wolf,  1980]: 

P  A  H  =  PL  +  PV 

The  weight  matrix  P13xl3,  coefficient  matrix  A13x9,  observation  matrix  L 
and  residual  matrix  V,3xl  are: 


'100  ^  f  1  0  0  0  0  0  0  0  0 

!  -1  10000000 

,  0-1  1000000 

1  ZEROS  00-1  100000 

1  0  0  0-1  1  0  0  0  0 

J  0000-1  1000 

,  00000-1  100 

!  A=  000000-1  10 
,  1  0  0  0  0  0  0-1  0 

,  001000-1  00 

ZEROS  1  0010000-1  0 

1  0  0  0-1  0  0  0  0  1 

,  00001000-1 

l  J  000-1  10000 

v  y 


The  program  Lobs. basic  [Jeyapalan,  1988]  was  used  to  find  H  for  the  temporary 
bench  marks.  The  program  (Appendix  B)  gives  the  H  and  V  matrices: 
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/  0.00000  \ 

0.00007 

/  0.000  \ 

0.00007 

-0.808 

-0.00022 

4.579 

-0.00044 

7.728 

-0.00022 

9.262 

V  = 

-0.00022 

8.246 

0.00007 

4.032 

-0.00019 

4.521 

-0.00010 

V  8.445  7 

-0.00010 

0.00022 

^0.00023  ) 

The  standard  deviation  of  the  observation  equations  is  0.36  mm. 

B.  GPS  DATA  PROCESSING 
1.  Automatic  Processing  Mode 

Trimble  recommends  that  an  automatic  processing  mode  be  used  if  more 
than  one  hour  of  data  from  four  or  five  GPS  satellites  is  taken  at  each  site  for 
baselines  lengths  up  to  30  km. 

The  output  files  contain  triple  difference,  double  difference  float  and 
double  difference  fixed  solutions.  The  double  difference  fixed  solution  provides 
the  most  accurate  solution  if  the  following  conditions  are  met: 

a.  Integer  bias  search  indicates  the  ratio  sum-of- squares  must  be  greater  than 
3.0; 

b.  RMS  (cycles)  less  than  0.05;  and 

c.  Difference  between  the  float  and  fixed  solutions  is  less  than  10  cm,  in  any 
component  of  the  baseline  (X,Y,or  Z). 

In  the  GPS  data  processing  all  the  solutions  from  the  bench  marks 
satisfied  these  conditions  except  station  J  697  for  which  RMS  was  0.068.  The 
distance  from  NPS  Building  224  to  J  697  is  33  km.  Trimble  recommends  a 
baseline  length  greater  than  30  km  when  the  fixed  solution  does  not  meet  the 
above  conditions.  The  float  solution  should  be  used  provided  the  RMS  of  fit  is 
better  than  0.08.  Ah  of  double  difference  fixed  solutions  were  used  for  the  bench 
marks  in  this  study.  The  double  difference  fixed  solutions  are  shown  in  Table  3. 


33 


I 


Table  3 

.  DOUBLE  E 

1IFFERE 

NCE  FIXED  SOLUTIONS 

Bench 

Slope  distance 

paga 

Coordinates  difference  between 

mark 

from  Bldg.  224 
(m) 

Ratio 

fixed  and  float  solution 

AX  (m)  AY  (m)  AZ  (m) 

TREE 

66 

I 

20.6 

-2.0 

-1.4 

GH1 

98 

257.3 

-0.5 

0.0 

GH2 

256 

0.025 

94.7 

-3.2 

0.3 

GH3 

536 

0.014 

276.4 

-1.1 

-0.1 

GH4 

547 

0.019 

151.0 

-0.7 

0.4 

-0.2 

GH5 

424 

0.014 

245.7 

0.0 

-1.6 

-1.2 

GH6 

166 

0.015 

254.9 

-0.7 

0.5 

-0.1 

GH7 

169 

0.015 

288.6 

-0.7 

0.3 

-0.1 

GH8 

413 

0.016 

160.0 

-0.4 

0.5 

0.1 

K  152 

14914 

0.035 

9.7 

-0.3 

0.1 

-0.4 

B  21 

2082 

0.014 

193.5 

0.4 

-0.2 

0.1 

S  812 

17062 

0.020 

15.8 

-1.2 

1.7 

0.1 

P  812 

20011 

0.020 

31.7 

0.9 

-1.4 

0.4 

D  697 

25249 

0.045 

11.1 

8.2 

9.3 

0.4 

J  697 

32692 

0.068 

6.3 

13.8 

-10.7 

2.6 

F  813 

16861 

0.022 

35.9 

7.1 

0.1 

GWM27 

16687 

0.049 

15.3 

4.4 

0.5 

Default  parameters  are  assumed  for  automatic  processing.  These 
parameters  have  a  minimum  elevation  mask  of  15°  for  double  differences. 
Station  one  coordinates  were  the  best  code  positions  during  data  logging.  Actual 
surface  meteorological  values  were  used  in  the  processing. 

For  sets  of  data  covering  more  than  one  hour,  epoch  increments  of  five 
or  ten  provide  full  precision.  Automatic  processing  begins  with  several  iterations 
of  triple  differences  to  establish  a  starting  value  for  the  double  difference 
solution.  A  cycle  slip  fixer  processes  every  epoch. 

Double  difference  solutions  begin  with  several  iterations  of  the  float 
solution  where  biases,  as  well  as  baseline  components,  are  computed.  It  is  called 
the  float  solution,  since  the  biases  are  allowed  to  float  and  are  not  constrained  to 
be  integers.  Each  iteration  should  indicate  convergence  toward  zero  for 
observed  minus  computed  phases.  The  default  elevation  mask  is  15°  for 
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processing  double  differences.  After  completion  of  the  float  solution, 
correlations  and  biases  are  computed.  The  integer  search  algorithm  begins  to 
determine  the  proper  integer  values  for  the  bias.  After  the  biases  are  set  to  their 
integer  values,  the  processing  is  repeated  to  compute  the  fixed  solution.  Ellipsoid 
height  differences  Ah  for  the  temporary  bench  marks  with  respect  to  NPS 
Building  224,  found  using  double  difference  fixed  solutions  is  given  in  Table  4. 


Table  4.  Ah  AND  WGS  84  COORDINATES  FOR  THE 
TEMPORARY  BENCH  MARKS. 


Bench  mark 

B1 

nnsm 

TREE 

-8.3532 

-2707298.015 

-4353450.286 

3781781.432 

GH1 

-9.1814 

-2707385.942 

-4353407.659 

3781790.139 

GH2 

-3.7766 

-2707502.139 

-4353548.202 

3781549.491 

GH3 

-0.6382 

-2707528.391 

-4353750.407 

3781314.047 

GH4 

0.9055 

-2707367.536 

-4353813.483 

3781306.605 

GH5 

-0.1099 

-2707202.581 

-4353799.531 

3781493.061 

GH6 

-4.3397 

-2707252.929 

-4353587.524 

3781666.956 

GH7 

-3.8509 

-2707406.088 

-4353553.901 

3781589.428 

GH8 

0.0843 

-2707329.452 

-4353759.387 

3781428.906 

Ah  for  the  permanent  bench  marks  with  respect  to  NPS  Building  224 
using  double  difference  fixed  solutions  are  given  in  Table  5. 
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Table  5.  Ah  AND  WGS  84  COORDINATES 
PERMANENT  BENCH  MARKS 


FOR  THE 


Bench  mark 

MESBMl 

\mBSsm 

wmmm 

K  152 

2.2398 

-2696890.306 

-4350931.132 

3792065.393 

B21 

-8.4275 

-2708436.436 

-4351973.362 

3782674.213 

S  812 

1.6335 

-2692137.111 

-4360870.528 

3784094.277 

P  812 

-1.3569 

-2689245.909 

-4360679.825 

3786312.651 

D  697 

15.1728 

-2685153.665 

-4357222.820 

3793175.171 

J  697 

71.3023 

-2679276.491 

-4356508.580 

3798216.777 

F  813 

-6.7070 

-2695613.779 

-4350454.843 

3793463.222 

GWM  27 

49.1051 

-2692408.721 

-4360427.654 

3784433.846 

Geoid  heights  of  bench  marks  were  obtained  by  using  the  Geoid.exe 
program  from  NGS  or  Trimvec  (Table  6). 


Table  6.  GEOID  HEIGHTS 


Bench  mark 


TREE  *34.7 

GH1  -34.7 

GH2  -34.7 

GH3  -34.7 

GH4  -34.7 

GH5  -34.7 

GH6  -34.7 

GH7  -34.7 

GH8  -34.7 

K  152  -34.1 

B  21  -34.7 

S  812  -33.9 

P  812  -33.8 

D  697  -33.6 

J  697  -33.3 

F  813  -34.0 

GWM  27  -33.9 


Trimvec  (m 


-34.30428 

-34.30671 

-34.31508 

-34.32072 

-34.31605 

-34.30745 

-34.30532 

-34.31140 

-34.31247 

-33.83303 

-34.31918 

-33.86473 

-33.76630 

-33.58265 

-33.42047 

-33.77899 

-33.86547 
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The  method  of  collocation  was  used  to  solve  the  local  geoid  model.  The 
procedure  is  described  as  follows  : 

a.  Reference  Station  Selection 

Station  TREE  was  selected  as  a  reference  station  on  the  NPS  campus. 
Coordinate  differences  AX,  AY  and  AZ  in  WGS  84,  for  the  temporary  bench 
mark  with  respect  to  the  TREE  were  computed.  These  are  shown  in  Table  7. 


Table  7.  AX,  AY  AND  AZ  ON  THE  NPS  CAMPUS 


Bench  mark 

wmm ■ 

AY  (m) 

wmm 

TREE 

0.000 

0.000 

0.000 

GH1 

-87.927 

42.627 

8.707 

GH2 

-204.124 

-97.916 

-231.941 

GH3 

-230.376 

-300.121 

-467.385 

GH4 

-69.521 

-363.197 

-474.827 

GH5 

95.434 

-349.245 

-288.371 

GH6 

45.085 

-137.238 

-114.476 

GH7 

-108.073 

-103.615 

-192.004 

GH8 

-31.437 

-309.101 

-352.526 

Station  K  152  was  selected  as  a  reference  station  in  the  large  off 
campus  area.  The  coordinate  differences,  AX,  AY  and  AZ  between  the 
permanent  bench  marks  and  K  152  were  computed  in  WGS  84  (Table  8). 
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Table  8.  AX,  AY  AND  AZ  IN  THE  OFF  CAMPUS  AREA 


Bench  mark 

AX  (m) 

■B9HSIH 

■■gam 

K  152 

0.000 

0.000 

0.000 

B  21 

-11546.130 

-1042.230 

-9391.180 

S  812 

4753.195 

-9939.396 

-7971.116 

P  812 

7644.397 

-9748.693 

-5752.742 

D  697 

11736.641 

-6291.688 

1109.778 

J  697 

17613.815 

-5577.448 

6151.384 

F813 

1276.527 

476.289 

1397.829 

GWM27 

4481.585 

-9496,522 

-7631.547 

b.  Determine  the  Local  Geoid  Model 

I  start  with  assuming  the  signal,  S,  equals  zero  and  the  weight  matrix 
P,  equals  the  unit  matrix.  Seven  control  marks  in  both  areas  led  to  seven 
observation  equations.  The  coordinate  differences  between  the  reference  marks, 
which  are  TREE  on  the  NPS  campus  and  K  152  in  the  off  campus  area,  are  the 
coefficients  of  the  parameters.  The  observation  equation  is  given  by 

AN  =  H  +  N0(X,Y,Z)  -  Ah 


where 

H  :  the  orthometric  height,  from  levelling 

N0(X,Y,Z)  :  the  global  geoid  height,  obtained  form  NGS  or  Trimvec 
Ah  :  the  ellipsoid  height  difference  with  respect  to  the  reference  station 
The  local  variation  of  the  geoid  height,  AN,  is  also  given  by 

AN  =  hQ  +  combination  of  coordinate  differences 

h0,  the  ellipsoid  height,  includes  a  constant  correction  to  N0(X,Y,Z)  at  the 
reference  station.  AN  for  temporary  bench  marks,  N0(X,Y,Z)  from  the  NGS 
(designated  as  NAN)  and  N0(X,Y,Z)  from  the  Trimvec  (designated  as  TAN)  are 
given  in  Table  9. 
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Table  9.  AN  ON  THE  NPS  CAMPUS 


Bench  mark 

TAN  (m) 

TREE 

-26.3468 

-25.9511 

GH1 

-26.3266 

-25.9333 

GH2 

-26.3444 

-25.9595 

GH3 

-26.3338 

-25.9545 

GH4 

-26.3435 

-25.9595 

GH5 

-26.3441 

-25.9515 

GH6 

-26.3283 

-25.9336 

AN  for  permanent  bench  marks  with  N0(X,Y,Z)  from  the  NGS  (designated  as 
NAN)  and  N0(X,Y,Z)  from  the  Trimvec  (designated  as  TAN)  are  given  in  Table 
10. 


Table  10.  AN  IN  THE  OFF  CAMPUS  AREA 


Bench  mark 

1 

TAN  (m) 

K  152 

-19.3748 

-19.1078 

B  21 

-20.4855 

-20.1047 

S  812 

-19.5555 

-19.5202 

P  812 

-19.3091 

-19.2754 

D  697 

-18.7158 

-18.6984 

J  697 

-19.5413 

-19.6618 

F  813 

-19.3100 

-19.0890 

A  least  squares  adjustment  program  Lobs. basic  was  used  to  determine 
the  geoid  model  which  has  the  smallest  standard  deviation.  This  was  done  by 
using  the  different  combinations  of  the  coordinate  differences.  GPS  data  from 
the  off  campus  area  were  studied  to  determine  the  local  geoid  model.  N0(X,Y,Z) 

from  the  Trimvec  program  were  used  to  compute  AN. 
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Four-parameter  combinations  and  their  standard  deviations  for  the  off 
campus  area  are  shown  in  Table  1 1 . 


Table  11.  FOUR  PARAMETERS 


Parameters 

Standard  deviation  (m) 

hQ  +  aAX  +  bAY  +  cAZ 

0.45 

ho  +  aAX  +  bAY  +  cAY2 

0.52 

hfl  +  aAX  +  bAY  +  cAX2 

0.34 

h0  +  aAX  +  bAY  +  cAXAY 

0.47 

h0  +  aAX  +  bAZ  +  cAY2 

0.52 

hfi  +  aAX  +  bAZ  +  cAX2 

0.34 

hfl  +  aAX  +  bAZ  +  cAXAY 

0.47 

hQ  +  aAX  +  bAY2  +  cAX2 

0.31 

h0  +  aAX  +  bAX2  +  cAXAY 

0.37 

hQ  +  aAY  +  bAZ  +  cAY2 

0.52 

h0  +  aAY  +  bAZ  +  cAX2 

0.34 

h0  +  aAY  +  bAZ  +  cAXAY 

0.46 

h0  +  aAY  +  bAY2  +  cAX2 

0.43 

hfl  +  aAY+ bAY2  +  cAXAY 

0.58 

h0  +  aAY+  bAX2  +  cAXAY 

0.28 

hQ  +  aAZ  +  bAY2  +  cAX2 

0.37 

hQ  +  aAZ  +  bAX2  +  cAXAY 

0.32 

hQ  +  aAZ  +  bAY2  +  cAXAY 

0.49 

h0  +  aAY2  +  bAX2  +  cAXAY 

0.20  * 

*  The  smallest  standard  deviation. 
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Five-parameter  combinations  and  their  standard  deviations  for  the  off  campus 
area  are  given  in  Table  12. 


Standard  deviation  (m 


0.43 


Table  12.  FIVE  PARAMETERS 


Parameters 


h  +  aAX  +  bAY  +  cAZ  +  dAX‘ 


h0  +  aAX  +  bAY  +  cAZ  +  dAY2 


h„  +  aAX  +  bAY  +  cAZ  +  dAXAY 


h0  +  aAX  +  bAZ  +  cAY2  +  dAXAY 


h  +  aAX  +  bAX2  +  cAY2  +  dAXAY 


h0  +  aAY  +  bAZ  +  cAX2  +  dAXAY 


h  +  aAY  +  bAZ  +  cAX2  +  dAXAY 


h  +  aAY  +  bAX2  +  cAY2  +  dAXAY 


h0  +  aAZ  +  bAX2  +  cAY2  +  dAXAY 


*  The  smallest  standard  deviation. 


Six-parameter  combinations  and  their  standard  deviations  for  the  off  campus  area 
are  given  in  Table  13. 


Table  13.  SIX  PARAMETERS 


Parameters 


hQ  +  aAX  +  bAY  +  cAZ  +  dAX2  +  eAXAY 


hQ  +  aAX  +  bAY  +  cAX2  +  dAY2  +  eAXAY 


h0  +  aAX  +  bAZ  +  cAX2  +  dAY2  +  eAXAY 


hQ  +  aAY  +  bAZ  +  cAX2  +  dAY2  +  eAXAY 


*  The  smallest  standard  deviation. 


Standard  deviation  (m 


0.26 
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c.  Evaluating  the  Geoid  Model  Using  Check  Marks 

In  order  to  evaluate  the  accuracy  of  the  geoid  model,  the  geoid  height 
model  can  be  rearranged  to  calculate  the  orthometric  height  of  check  marks 
which  are  GH7  and  GH8  on  the  NPS  campus  and  GWM  27  in  the  off  campus 
area. 

H  =  h0+  Ah  -  N0(X,Y,Z)  +  aAY  +  bAX2  +  cAY2  +  dAXAY 
or 

H  =  h0+  Ah  -  N0(X,Y,Z)  +  a’AX  +  b’AZ  +  c’AX2  +  d’AY2  +  e’AXAY 

The  h0,  a,  b,  c  and  d  of  the  five-parameter  geoid  model,  and  h0,  a',  b',  c',  d’  and 

e'  of  the  six -parameter  geoid  model,  can  be  solved  by  least  squares  adjustments. 
The  results  for  h0,  a,  b,  c  and  d  for  the  five-parameter  model  using  the  seven 
control  marks  in  the  off  campus  area  are  given  in  Table  15.  The  results  for  h0, 
a,  b,  c  and  d  are  shown  in  the  NGS  column,  where  the  N0(X,Y,Z)  was  from  the 
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NGS  Geoid.exe  program.  The  results  for  h0,  a,  b,  c,  d  and  standard  deviation,  a, 
are  shown  in  the  Trimvec  column,  where  the  N0(X,Y,Z)  was  from  the  Trimvec 
Geoid.exe  program. 


Table  15.  h0,  a,  b,  c,  d,  ct  FOR  SEVEN  CONTROL  MARKS 


Parameter 


ho 


-19.26895 
-2.344959E-04 
-8.523017E-09 
-3.5743 14E-08 
-2.211550E-08 
0.130518 


Trimvec  (m 


-19.00915 

-2.736598E-04 

-8.528133E-09 

-3.905052E-08 

-1.728313E-08 

0.095723 


The  results  of  h0,  a,  b,  c,  d  and  o  of  the  five-parameter  geoid  model 
using  the  six  control  marks  (excluding  B  21)  in  the  off  campus  area  are  given  in 
Table  16. 


Table  16.  h0,  a,  b,  c,  d,  o  FOR  SIX  CONTROL  MARKS 


■19.25415 

-2.752543E-04 

-7.735253E-09 

-3.777950E-08 

-1.798617E-08 

0.171919 


Trimvec  (m 


-19.01411 

-2.600850E-04 

-8.792540E-09 

-3.836067E-08 

-1.867193E-08 

0.133516 


The  results  of  h0,  a,  b,  c,  d  and  a  of  the  five-parameter  geoid  model 
using  the  six  control  marks  (excluding  J  697)  in  the  off  campus  area  are  given  in 
Table  17. 


} 


Table  17.  h0,  a,  b,  c,  d,  q  FOR  SIX  CONTROL  MARKS 


Parameter 

NGS  (m) 

Trimvec  (m) 

ho 

-19.265355 

-19.03487 

a 

-2.503395E-04 

-1.977533E-04 

b 

-8.731945E-09 

-7.5264 10E-09 

c 

-3.711466E-08 

-3.246532E-08 

d 

-2.179058E-08 

-1. 88279  IE-08 

CT 

0.1841317 

0.1206984 

The  results  of  1^,  a,  b,  c  and  d  of  the  five-parameter  geoid  model  using 
the  five  control  marks  (excluding  B  21  and  J  697)  in  the  off  campus  area  are 
given  in  Table  18.  The  standard  deviation  is  undefined  since  degree  of  freedom 
equals  zero. 


Table  18.  h0,  a,  b,  c,  d  FOR  FIVE  CONTROL  MARKS 


Parameter 

NGS  (m) 

Trimvec  (m) 

ho 

-19.37497 

-19.10791 

a 

1.032352E-04 

3.397465E-05 

b 

8.119969E-09 

3.521564E-09 

c 

7.37782  IE-09 

-3.339665E-09 

d 

1.338776E-09 

-3.667083E-09 

The  results  of  h0,  a,  b,  c,  d  and  ct  of  the  five-parameter  geoid  model 
using  the  seven  control  marks  on  the  NPS  campus  are  given  in  Table  19. 


Table  19.  h0,  a,  b,  c,  d,  ct  FOR  SEVEN  CONTROL  MARKS 


Parameter 

NGS  (m) 

Trimvec  (m) 

ho 

-26.33543 

-25.94105 

a 

2.976507E-06 

1.192093E-07 

b 

-5.419133E-08 

-1.763692E-07 

c 

-3.601599E-08 

-8.489588E-08 

d 

6.30971  IE-08 

-3.702007E-08 

CT 

0.013773 

0.014612 
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The  results  of  l^,  a',  b',  c',  d',  e'  and  a  of  the  six-parameter  geoid  model 
using  the  seven  control  marks  in  the  off  campus  area  are  shown  in  Table  20. 

Table  20.  hfl,  a’,  b’,  c',  d’,  e',  o  FOR  SEVEN  CONTROL  MARKS 


Parameter 


h« 


Trimvec  (m 


-19.27388 

-19.01660 

2.083 182E-04 

1.842380E-04 

-2.561510E-04 

-2.520234E-04 

-7.731387E-09 

-8.335974E-09 

-3.767036E-08 

-3.960304E-08 

-1.249646E-08 

-1.521767E-08 

0.137149 

0.123960 

The  results  of  h0,  a',  b',  c',  d'  and  e'  of  the  six-parameter  geoid  model 
using  the  six  control  marks  (excluding  B  21)  in  the  off  campus  area  are  given  in 
Table  21.  The  standard  deviation  is  undefined  since  degree  of  freedom  equals 
zero. 

Table  21.  hfl,  a',  b',  c',  d',  e'  FOR  SIX  CONTROL  MARKS 


Parameter 


hn 


-19.37489 
2.57641  IE-04 
-1.691282E-04 
-1.083072E-08 
-2.817797E-08 
-5.824404E-09 


Trimvec  (m 


-19.10785 
2.288222E-04 
-1.734756E-04 
-1.113654E-08 
-3.1 0274  IE-08 
-9.185896E-09 


The  results  of  h0,  a',  b',  c\  d'  and  e'  of  the  six-parameter  geoid  model 
using  the  six  control  marks  (excluding  J  697)  in  the  off  campus  area  are  given  in 
Table  22.  The  standard  deviation  is  undefined  since  degree  of  freedom  equals 
zero. 
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Table  22.  h0,  a*,  b’,  c’,  d*,  e*  FOR  SIX  CONTROL  MARKS 


Parameter 

NGS  (m) 

Trimvec  (m) 

ho 

-19.37466 

-19.10774 

a’ 

1.065433E-04 

9.238720E-05 

b' 

-4.604459E-05 

-6.231666E-05 

c’ 

-2.0127 12E-09 

-3.170499E-09 

d’ 

-1.141234E-08 

-1.586523E-08 

e’ 

-2.528395E-09 

-6.206392E-09 

The  results  of  h^,,  a’,  b’,  c',  d',  e'  and  a  of  the  six-parameter  geoid  model 
using  the  seven  control  marks  on  the  NPS  campus  are  given  in  Table  23. 

Table  23.  h0,  a',  b',  c*,  d',  e\  a  FOR  SEVEN  CONTROL  MARKS 


Parameter 

NGS  (m) 

Trimvec  (m) 

ho 

-26.33389 

-25.93820 

a’ 

1.645088E-05 

4.556775E-05 

b' 

4.225970E-05 

6.169081E-05 

c* 

6.856863E-08 

6.391201E-08 

d’ 

6.00703  IE-08 

5.878974E-08 

e' 

1.781 154E-07 

1.767185E-07 

o 

0.018859 

0.018872 

Fortran  program  GPSDIS  (Appendix  C)  was  used  to  calculate  the 
orthometric  heights  of  the  check  marks.  The  results  are  discussed  in  the  Chapter 
VII. 

d.  Obtaining  the  Weight  Matrix 

Since  the  standard  error  of  the  geoid  model  is  about  ±  2  cm,  which  is 
the  same  or  better  than  GPS  Ah  accuracy,  it  is  not  actually  necessary  to  do 
further  computation.  Nevertheless,  the  weight  matrices,  P,  were  computed. 

The  five-parameter  geoid  model  used  for  the  NPS  campus  and  the 
six-parameter  geoid  model  used  for  the  off  campus  area  were  employed  in 
obtaining  P.  The  covariance  function,  C(r),  of  the  control  marks  (p.  18)  can  be 
obtained  by  using  the  least  squares  residuals.  The  coefficients  B  and  k  of  C(r) 
can  be  obtained  by  solving  simultaneously  the  equations  for  C(r)  with  different 
distances  from  the  centers. 


(1)  P  for  NPS  Campus.  The  residuals  of  the  five-parameter  geoid 
model  of  the  control  marks  on  the  NPS  campus  are  given  in  Table  24.  The 
global  geoid  heights  were  from  the  NGS. 


Table  24.  RESIDUALS  FOR  FIVE-PARAMETER  MODEL 


Control  mark 

#  of  residual 

Residual  (m) 

TREE 

v, 

1.137352E-02 

GH1 

v2 

-9.420395E-03 

GH2 

V3 

7.339478E-03 

GH3 

V4 

-4.278 183E-03 

GH4 

V5 

3.572464E-03 

GH5 

V6 

6.446839E-04 

GH6 

V7 

-8.71 467 6E-03 

For  Tj  =  593  m 


C(ri)  = 


v  v  +  V  V  +  V  V  +  v  V  +  V  V, 

14  15  24  25  46 


=  -2.194272E-03 


For  r2  =  319  m 


V1V3+V2VV2V7+V,VV,VS+V.Vfi+V.V7+V4V,+VSV, 

C(r  )=  — L- 2 — 2-2 — 2_Z — 2_i — — 2— £ — 2—1 — 4_$ — S_&  =  3 .945813E-6 
2  5 


Thus  the  equations  of  covariance  can  be  written  as 

C(rj)  =0.0138 +  B  sinCkTj) 
and 

C(r2)  =  0.0138  +  B  sin(kr2) 

Approximate  solutions  for  coefficients  B  and  k  can  be  found  by  using  the 
Maclaurin  series  expansion.  The  solved  covariance  function  is  given  by 

C(r)  =  0.0138  -  0.0160  sin(1.7670  r) 
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Cq  can  be  obtained  by  using  the  C(r)  with  the  distances  between  the  control 

marks.  Fortran  program  DISTCO  (Appendix  D)  was  used  to  compute  the 

distances  between  the  control  marks,  as  well  as  Cq  and,  hence,  the  weight  matrix,  | 

P  =  C^1 .  The  Matlab  program  on  the  NPS  mainframe  computer  was  used  to  find 

the  inverse  of  Cq  (P  =  Cq1).  P  for  the  five-parameter  geoid  model  for  the  NPS 

campus  is  given  by 


P  = 


\ 


162 


ZEROS 


142 


78 


98 


102 


ZEROS 


87 


100 


(2)  P  for  off  Campus  Area.  The  residuals  of  the  six-parameter  geoid 
model  of  the  control  marks  in  the  off  campus  area  are  given  in  Table  25.  The 
global  geoid  heights  were  calculated  using  the  Trimvec  program. 


Table  25.  RESIDUALS  FOR  SIX-PARAMETER  MODEL 


Control  mark 

#  of  residual 

Residual  (m) 

K  152 

v. 

9.120178E-02 

B  21 

V2 

-9.773254E-03 

S  812 

V3 

6.391526E-03 

P  812 

V4 

1.983643E-04 

D  697 

V5 

-2.77996  IE-02 

J  697 

V6 

1.685524E-02 

F  813 

v. 

-7.65171  IE-02 

For  r  =  13994  m 


C(fl)  = 


V.V  +  V,V7  +  V,V4  +  V,V<  +  V7V7  +  V.V7 
- 2—2 - L_3 - 2-  4-  15 - 3_7 - 4__7_  _  -6.659883E-04 
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For  r2  =  19230  m 


C(r2)  = 


v  v  +V  V  +V  V 

1  6  2  3^  v3  v6 

2 


7.912463E-04 


Thus  the  equations  of  covariance  can  be  written  as 

CO-j)  =  0.0124  +  B  sin^) 
and 

C(r2)  =  0.0124  +  B  sin(kr2) 

Approximate  solutions  for  coefficients  B  and  k  can  be  found  by  using  a 
Maclaurin  series  expansion.  The  solved  covariance  function  is  given  by 

C(r)  =  0.0124  -  0.1329  sin(2.8278  r) 

Similarly,  Cq  can  be  obtained  by  using  the  C(r)  with  the  distances  between  the 

control  marks.  Fortran  program  DISTCO  (Appendix  D)  was  used  to  compute 
the  distances  between  the  control  marks  and  Cq,  to  obtain  the  weight  matrix  (P  = 
Cq'1).  The  Matlab  program  on  the  NPS  mainframe  computer  was  used  to  find  the 
inverse  of  Cq  (P  =  C^1).  P  for  the  six-parameter  geoid  model  of  the  off  campus 
area  is  given  by 


r 


p  = 


25 


15 

15 

ZEROS 


ZEROS 


16 


10 


17 


25 


V 
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VI. 


Factors  affecting  the  GPS  measured  ellipsoid  height  differences,  Ah, 
include:  (1)  comparing  GPS  differencing  solutions,  (2)  standard  error  of  GPS 
observations,  (3)  corrections  for  surface  meteorological  values,  and  (4) 
observation  durations  for  GPS. 

A.  COMPARING  GPS  DIFFERENCING  SOLUTIONS 

1.  Ellipsoid  Height  Differences  Tested  at  a  Fixed  Position 

In  order  to  determine  the  best  solution  for  Ah  using  GPS,  the  Ah  data 
from  a  Trimble  4000SX  GPS  receiver  were  collected  in  the  K  parking  lot  on  the 
NPS  campus  from  0633  to  1007  UTC  on  February  19,  1988.  The  design  of  the 
experiments  is  shown  in  the  Figure  10. 


Figure  10.  Illustration  of  the  Ah  test  of  Trimble  4000SX  GPS 


receiver. 
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A  plumb  bob  was  used  to  point  to  a  fixed  position  (36°  35'  56"  N,  121° 
52’  36”  W)  on  the  ground.  The  observations  were  taken  by  changing  the  antenna 
height  over  a  distance  of  about  60  cm,  which  was  measured  with  a  tape  every 
half  hour.  The  ellipsoid  height  differences.  Ah,  and  the  orthometric  height 
differences,  AH,  should  agree  with  the  change  of  antenna  height  for  a  fixed 
position.  Triple  difference  solutions,  Ahyjy,  are  given  in  Table  26.  Table  also 
gives  the  differences  between  the  Ah^  at  succeeding  observation  times,  DAhTRI; 
double  difference  float  solutions,  AhFLT;  the  differences  between  the  Ah^  at 
succeeding  observation  times,  DAhFLT;  and  double  difference  fixed  solutions, 
AhFK,  and  the  differences  between  the  Ah^  at  succeeding  observation  times, 
DAh  Fix- 


Table  26.  Ah  AND  DAh  FOR  DIFFERENT  ANTENNA  HEIGHTS 


A^tri 

DAhjxj 

AhpLT 

DAhpLT 

Ah  fix 

DAhpjx 

(19-2-1988) 

(m) 

(m) 

(m) 

(m) 

(m) 

(m) 

0633  -  0712 

4.4410 

-0.7189 

4.4429 

-0.6339 

4.5045 

-0.5463 

0736  -  0808 

5.1599 

-0.3973 

5.0768 

-0.5676 

5.0508 

-0.5964 

0824  -  0855 

5.5572 

-0.1954 

5.6444 

-0.7144 

5.6472 

-0.7030 

0911  -0943 

5.7526 

-0.4531 

6.3588 

-0.2920 

6.3502 

-0.2660 

0945  -  1007 

6.2057 

6.6508 

6.6162 

The  differences  between  the  GPS  DAh  and  the  differences  of  the 
measured  ellipsoid  height  differences,  DAhMEA,  are  given  in  Table  27. 


51 


Table  27.  THE  DIFFERENCE  OF  Ah 


ANTENNA 
HEIGHT  (m) 

Eggai 

fMGSM 

■HSHi 

■HSM 

2.809 

2.262 

-0.5470 

0.1719 

0.0869 

-0.0007 

1.666 

-0.5960 

-0.1987 

-0.0284 

0.0004 

1.074 

-0.5920 

-0.3966 

0.1224 

0.1110 

0.484 

-0.5900 

-0.1369 

-0.2980 

-0.3240 

The  results  show  that  the  double  difference  fixed  solution  is  the  best. 
The  difference  between  the  measured  ellipsoid  height  difference  and  the  GPS 
ellipsoid  height  difference  was  about  1  mm.  The  last  two  results  were  not  good 
because  they  had  multipath  effects  and  imaging  effects  from  the  ground  (antenna 
heights  1.074  m  and  0.484  m),  so  the  solutions  of  AH  compared  to  the  true 
differences  were  large.  This  indicates  that  the  antenna  should  be  at  least  1  m 
above  the  ground. 

2.  Ellipsoid  Height  Differences  Tested  at  Different  Positions 

Since  the  geoid  on  the  NPS  campus  has  a  flat  geoid  slope,  the  differences 
between  Ah,  DAh,  and  the  orthometric  height  difference,  AH,  of  the  temporary 
bench  marks  can  be  neglected.  The  Ah  and  DAh  from  the  GPS  observations  are 
given  in  Table  28.  The  differences  between  the  GPS  DAh  and  the  AH  are  given 
in  Table  29. 
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Table  28. 


Bench 

mark 


Ah  AND  DAh  ON  THE  NPS  CAMPUS 


8.3358 

-9.1859 

-3.8387 

-0.6393 


0.8848 


-0.0929 

-4.3584 


-3.8652 


0.0760 


0.8501 


-5.3472 


-3.1994 


-1.5241 

0.9777 

4.2655 

-0.4932 


-3.9412 


-8.3375 

-9.1688 

-3.7998 

-0.6374 


0.9067 


-0.1104 

-4.3407 


-3.8491 


0.0869 


0.8313 

-5.3690 


-3.1624 


-1.5441 


1.0171 

4.2303 

-0.4916 

-3.9360 


-8.3532 


-9.1814 

-3.7766 


-0.6382 


0.9055 

-0.1099 

-4.3397 

-3.8509 


0.8282 


-5.4048 

-3.1384 

-1.5437 

1.0154 


4.2298 


-0.4888 

-3.9352 


0.0843 


Table  29.  DIFFERENCES  BETWEEN  DAh  AND  AH 


Bench 


mark 


TREE  0.000 

GH1  -0.808 

GH2  4.579 

GH3  7.728 


9.262 


8.246 


4.032 


4.521 

8.445 


0.808 


-5.387 


-3.149 


-1.534 

1.016 


4.214 


-0.489 


-3.924 


-0.0423 

-0.0398 

0.0504 

-0.0099 

0.0383 

-0.0515 

0.0042 

0.0172 


-0.0235 

-0.0180 

0.0134 

0.0101 


-0.0011 


-0.0163 


0.0026 


0.0120 


-0.0204 

0.0178 

-0.0106 

0.0097 


0.0006 


-0.0158 


-0.0002 


0.0112 


It  is  clear  that  the  double  difference  fixed  solutions  give  the  best  solution 
in  a  flat  geoid  slope  area. 

B.  STANDARD  ERROR  OF  GPS  OBSERVATIONS 

To  determine  die  accuracy  of  GPS,  measurements  were  taken  at  bench  mark 
TREE  on  the  NPS  campus.  Three  measurements,  each  of  about  three  hours 
duration,  were  taken  for  this  analysis.  Double  difference  fixed  solutions  were 
used  to  analyze  the  data.  The  double  difference  fixed  solutions  for  these 
measurements  are  given  in  Table  30. 


Table  30.  Ah  AT  STATION  TREE 


Date 

Time  (UTC) 

AhFIX  (m) 

Feb.  5,  88 

0730  -  1025 

-8.4262 

Feb.  7,  88 

0721  -  1004 

-8.3532 

Feb.  11,  88 

0705  -  0959 

-8.3602 

The  standard  error  of  the  mean  of  Ah  (a  /  Vn)  is  about  ±  2  cm  at  TREE  and 
the  standard  error  (o)  of  the  single  measurement  is  about  ±  4  cm  for  three  hours 
of  observation. 

Results  obtained  by  Strange  in  1985  for  a  project  southeast  of  Phoenix, 
Arizona,  show  that  using  multiple  GPS  observations  at  the  same  point  lead  to 
uncertainties  typically  less  than  2  cm  in  Z-direction  over  20-km  distances 
[Zilkoski,  1988]. 

Thus  the  accuracy  of  the  ellipsoid  differences  obtainable  by  GPS 
measurements  is  expected  to  be  about  ±  2  cm  for  about  three  hours  of 
observation. 

C.  CORRECTIONS  FOR  SURFACE  METEOROLOGICAL 
VALUES 

The  meteorological  correction  is  a  function  of  the  atmospheric  refractivity 
as  computed  by  the  surface  meteorological  values  of  pressure,  temperature, 
humidity  and  the  elevation  angle  of  the  satellite.  Larger  corrections  are  required 
for  low  elevation  angles,  as  the  signal  travels  a  longer  path  through  the 
troposphere.  A  modified  Hopfield  troposphere  model  [Fell,  1975]  is  used  for 
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r 

t 

i 

|  automatic  processing  of  Trim640  software.  The  model  corrects  for  at  least  90% 

of  the  tropospheric  delay  [Remondi,  1984].  Default  surface  meteorological 
pH  values  for  automatic  processing  are  1010  millibars,  20°  C  and  50  %  relative 

[  humidity.  These  parameters  can  be  reset  before  the  processing.  Trim640  allows 

|  only  one  pressure,  temperature  and  humidity  entry  for  each  site  per  session. 

I  Ellipsoid  height  differences,  Ahp^,  were  calculated  for  four  marks  on  the  NPS 

r  campus  and  two  off  campus  using  the  default  meteorological  parameters  and 

using  the  true  values.  There  are  given  in  Table  31  along  with  differences  found 
[  between  the  two  methods. 


Table  31.  COMPARISON  OF  SURFACE  METEOROLOGICAL 
CORRECTION 


Bench  mark 

Default  Ahpjx  (m) 

Difference  (mm) 

GH2 

-3.7777 

-3.7766 

1.1 

GH3 

-0.6376 

-0.6383 

0.6 

GH4 

0.9070 

0.9055 

1.5 

GH7 

-3.8524 

-3.8509 

1.5 

K  152 

2.2398 

2.2430 

-3.2 

J  697 

-71.3179 

-71.3023 

15.6 

The  difference  is  about  ±  2  mm  on  the  NPS  campus  and  about  ±  2  cm  in  the 
Monterey  Bay  area.  The  Trim640  program  contains  no  ionospheric  model.  For 
Pi  first-order  relative  positioning  (1  part  in  105)  and  a  baseline  shorter  than  50  km 

ionospheric  delay  can  be  neglected  [Remondi,  1984],  In  this  study  surface 
meteorological  corrections  were  made  during  data  processing  to  get  the  best  Ah 
solutions. 


D.  OBSERVATION  DURATIONS  FOR  GPS 

To  find  out  how  long  observations  must  be  made  to  obtain  the  best  solutions 
in  the  field  when  batteries  are  used  for  power,  the  GPS  data  from  S  812,  J  697 
and  K  152  were  segmented  into  averaging  periods  of  length  nAt,  where  At  =  10 
minutes  and  n  =  1,  2,  3,  ...,  10.  Default  meteorological  values  were  used  for  this 
processing.  Ah^  for  various  observation  durations  for  the  stations  is  given  in 

Table  32  and  plotted  in  Figure  1 1  to  Figure  13. 
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Table  32.  Ah_  AS  A  FUNCTION  OF  OBSERVATION  DURATION 


I^E3SSE!!!S9^HB333SflB 

wnsxsm 

10 

1.5309 

71.5291 

2.3135 

20 

1.6141 

70.9697 

2.3102 

30 

1.6441 

71.2185 

2.3062 

40 

1.6539 

71.4514 

2.3065 

50 

1.6558 

71.4762 

2.3060 

60 

1.6498 

71.5343 

2.2563 

70 

1.6464 

71.5309 

2.2504 

80 

1.6349 

71.3407 

2.2504 

90 

71.3179 

2.2548 

100 

2.2430 

Figure  11.  Ah  vs.  observation  durations  at  station  S  812. 
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Figure  12.  Ah  vs.  observation  durations  at  station  J  697. 


Figure  13.  Ah  vs.  observation  durations  at  station  K  152. 

It  is  clear  that  Ah  was  affected  by  the  observation  duration.  Ah  varied  quite 
a  lot  in  the  first  forty  minutes.  After  about  forty  minutes.  Ah  tended  to  close  to 
the  mean  value  of  the  observations.  After  sixty  minutes.  Ah  didn’t  vary  much. 
Trimble  recommends,  to  ensure  sufficient  quality  data  and  to  obtain  first-order 
results  on  single  observations,  at  least  one  hour  should  be  taken  for  a  four  or  five 
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satellite  session.  For  baseline  lengths  up  to  30  km  the  conditions  for  the  fixed 
double  solution  of  automatic  processing  will  be  met  if  data  are  collected  for 
ninety  minutes  on  four  or  five  satellites  [Trimble,  1987b].  The  results  shows  that 
the  best  solutions  can  be  improved  with  longer  observation  times.  Normally,  at 
least  sixty  minutes  of  observations  are  required  for  vertical  control. 


VII. 


A.  EVALUATATION  OF  GEOID  MODEL 

To  find  the  accuracy  of  the  geoid  model  the  orthometric  heights  of  the 
check  marks  from  GPS,  H^g,  were  calculated,  and  these  were  compared  to  the 
orthometric  height  of  the  check  marks  from  the  levelling,  H,  FVFI ,  ^ 

I  started  with  the  assumptions  that  Cq  =  I  and  CM  =  0.  Five-parameter  and 
six-parameter  geoid  models  were  used  to  calculate  the  orthometric  height  of  the 
check  marks.  GWM  27  is  a  check  mark  in  the  off  campus  area.  GH7  and  GH8 
are  the  check  marks  on  the  NPS  campus.  N0(X,Y,Z)  was  calculated  using  both 

the  NGS  and  the  Trimvec  programs. 

1.  from  Five-Parameter  Geoid  Model 
a.  Hqps  in  the  off  Campus  Area 

Hgps  in  GWM  27  was  calculated  by  using  seven  control  marks,  six 
control  marks  (excluding  B  21  or  J  697)  and  five  control  marks  (excluding  B  21 
and  J  697).  Comparisons  of  H  in  GWM  27  are  given  in  Table  33. 


Table  33.  COMPARISONS  OF  H  USING  FIVE-PARAMETER 
MODEL  AT  GWM  27 


GWM  27 

^GPS  '  ^LEVELLING  (cm) 

#  of  control  marks 

Trimvec 

NGS 

7 

-11.55 

-2.24 

6  (excluding  B21) 

-9.77 

-8.07 

6  (excluding  J  697) 

-4.87 

-3.65 

5 

2.14 

6.58 

b.  Hgpj  on  the  NPS  Campus 

Heps  in  GH7  and  GH8  were  calculated  by  using  seven  control  marks. 
Comparisons  of  H  in  GH7  and  GH8  are  given  in  Table  34. 
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Table  34.  COMPARISONS  OF  H  USING  FIVE-PARAMETER 
MODEL  AT  GH7  AND  GH8 


Five  parameters 


Check  mark  |  Trimvec 


^GPS  ~  ^LEVELLING  (cm) 


.52 

0.82 

.20 

0.00 

a.  HqPS  in  the  off  Campus  Area 

Hqps  in  GWM  27  were  calculated  by  using  seven  control  marks,  and 
six  control  marks  (excluding  B  21  or  J  697).  Comparisons  of  H  at  GWM  27  are 
given  in  Table  35. 

Table  35.  COMPARISONS  OF  H  USING  SIX-PARAMETER 
MODEL  AT  GWM  27 


GWM  27 


#  of  control  marks 


6  (excluding  B21) 


b.  Hqps  on  the  NPS  Campus 

Hqps  of  GH7  and  GH8  were  calculated  by  using  seven  control  marks. 
Comparisons  of  H  at  GH7  and  GH8  are  given  in  Table  36. 

Table  36.  COMPARISONS  OF  H  USING  SIX-PARAMETER 
MODEL  AT  GH7  AND  GH8 


Six  parameters _ Hgps  -  HLEVELLING  (cm) 


Check  mark  Trimvec  NGS 


GH7  1.22  1.22 

GH8  0.22  0.25 


^GPS  ‘  ^IJFVF.I 1  TNG  (cm) 

Trimvec 

NGS 

-12.40 

-11.13 

-9.37 

-7.72 

-3.49 

-1.01 

Since  the  accuracy  of  predicted  geoid  height  at  check  marks  (±  0.2  to  2 
cm)  is  less  than  the  accuracy  of  GPS  observations  (±  2  to  4  cm),  it  may  be 
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concluded  that  P  =  I  is  a  satisfactory  covariance  matrix  and  further  computation 
is  unnecessary. 

B.  ACCURACY 

The  best  accuracy  of  the  five-parameter  geoid  model  is  ±  2  cm  in  the  off 
campus  area  and  is  ±  2  to  10  mm  for  the  NPS  campus.  The  best  accuracy  of  six- 
parameter  geoid  model  is  ±  1  cm  in  the  off  campus  area  and  is  ±  2  to  10  mm  on 
the  NPS  campus. 

C.  DISCUSSION 

The  errors  in  the  geoid  model  include  errors  associated  with  GPS  and  with 
levelling. 

1.  Errors  Associated  With  GPS 

a.  Atmospheric  Delays 

The  time  delay  of  RF  signals  passing  through  the  ionosphere  is  due  to 
a  reduction  in  speed  and  the  bending  of  the  ray,  both  effects  being  due  to 
refraction.  The  overall  delay  in  the  signal  is  nearly  inversely  proportional  to  the 
square  of  the  frequency.  The  ionospheric  delay  may  be  calculated  by  using  two- 
frequency  receivers  (C/A  code,  P  code),  and  corrected  for  with  a  satisfactory 
degree  of  accuracy.  Tropospheric  delays  are  independent  of  frequency.  They 
are  relatively  small  and  can  be  modelled  using  the  elevation  angle  of  the  satellite. 
The  Trimble  4000SX  receiver  only  observes  the  C/A  code  .  The  Trim640 
program  contains  no  ionospheric  modelling.  The  modified  Hopfield  model  [Fell, 
1975]  is  used  for  the  tropospheric  delay.  The  combined  effects  of  unmodelled 
ionospheric  and  tropospheric  errors  are  estimated  to  result  in  a  satellite-to-user 
range  error  of  from  2.4  to  5.2  m  [Milliken,  1980]. 

b.  Selection  of  Appropriate  Control  Marks 

Only  eight  permanent  bench  marks  were  recovered  in  the  off  campus 
study  area.  D  697  was  obstructed  by  a  tree.  F  813  was  set  near  a  steel  bridge. 
Although  offset  levelling  was  performed,  the  ellipsoid  height  could  have  been 
affected  by  offsetting.  The  shape  of  the  bench  marks  were  spread  out  roughly  in 
a  rhomboid  shape.  B  21  and  J  697  are  on  the  ends  of  the  rhomboid.  The 
distance  is  about  33  km  from  B  21  to  J  697.  The  elevation  changes  form  5.787  to 
85.061  m.  If  more  bench  marks  could  be  recovered  and  be  connected  to  these 
marks,  the  accuracy  of  the  seven  control  marks  could  be  improved. 
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c.  Error  of  Antenna  Height  Measurements  j 

The  slant  heights  of  GPS  antenna  were  measured  by  using  a  steel  tape  J 

before  and  after  data  collection.  Thus,  the  heights  of  antenna  were  calculated.  *1 

The  systematic  error  of  the  tape,  reading  error,  calculation  error  and  a  slack  tape 
doing  measurement  will  affect  the  accuracy  of  the  ellipsoid  height  difference. 

The  error  of  antenna  height  measurement  is  ±  2  to  10  mm. 

d.  Selecting  observation  period  with  optimal  satellite  geometry  A 

Because  there  are  only  seven  GPS  satellites  in  operation  at  present, 

the  observing  window  is  only  four  to  five  hours  a  day  for  four  or  five  satellites. 

The  observation  schedule  should  be  planed  early  and  considered  the  travel  time, 

and  time  to  reoccupy  the  marks.  In  this  study  most  of  the  GPS  measurements  Jj 

were  made  at  night.  PDOP  were  not  always  optimum. 

e.  Error  in  the  Computing  of  Trim640  Program 

There  is  no  standard  specification  of  the  GPS  software  at  present. 

Round  off  error  in  the  data  processing,  the  data  record  frequency,  and  the 
precision  of  the  ephemeris  could  affect  the  components  of  the  baseline. 

Specifications  of  the  GPS  software  certainly  will  be  developed  in  the  future. 

2.  Errors  in  Levelling 

The  error  in  levelling  would  cause  the  errors  of  the  observation  « 

equations  when  determining  the  local  geoid  model.  The  maximum  allowable 
closing  error  between  the  forward  and  backward  running  of  a  section  is  3.0  mm 
Vk  for  first-order  class  I  levelling,  and  9.0  mm  Vk  for  third-order  class  I 
levelling,  where  k  is  the  length  of  the  section  measured  in  km  [Bodnar,  1975].  ■ 

The  average  distance  from  K  152  to  the  permanent  bench  marks  is  about  13  km 
in  the  off  campus  study  area.  Thus,  for  first-order  class  I  levelling,  the 
maximum  allowable  closing  error  would  be  about  1 1  mm  in  the  off  campus  area. 

The  first-order  bench  marks  were  set  around  1933.  The  elevations  of  the  marks  J 

may  have  changed  caused  by  subsidences  or  earthquakes.  Since  the  time  did  not 
permit  rerunning  the  levelling  for  the  off  campus  area,  published  elevations  were 
assumed  correct.  Errors  in  levelling  could  contribute  to  errors  in  the  geoid 
model. 

The  distance  of  the  interlocking  levelling  circuit  on  the  NPS  campus  is 
about  0.8  km.  For  third-order  class  I  levelling,  the  maximum  allowable  closing 
error  is  about  7  mm  on  the  NPS  campus.  The  5  mm  closing  error  was  obtained 
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for  that  region.  The  geoid  height  model  thus  had  a  better  accuracy  on  the  NPS 
campus. 
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VIII. 


Miacnxviicija 


Local  geoid  models  for  the  Monterey,  California,  area  were  determined  by 
GPS  differential  positioning  using  the  method  of  collocation.  The  local  geoid 
models  were  based  on  Rapp's  360  degree  x  360  order  global  geoid  models 
determined  from  gravity  measurements.  Control  data  were  adjusted  by  least 
squares  to  solve  for  the  parameters  in  the  local  geoid  model.  Also  studied  were 
factors  that  affected  the  GPS-measured  ellipsoid  height  differences.  These  included 
(1)  comparing  GPS  differencing  solutions,  (2)  standard  error  of  GPS  observations, 
(3)  corrections  for  surface  meteorological  values,  and  (4)  observation  durations  for 
GPS.  Conclusions  are  the  following: 

1 .  The  accuracy  of  a  local  geoid  height,  N,  is  the  same,  whether  the  NGS  or 
the  Trimvec  version  of  Rapp's  global  geoid  model  is  used,  although  N0(X,Y,Z) 

calculated  using  Trimble's  Trimvec  program  gives  the  best  results  for  points  on  the 
NPS  campus,  and  N0(X,Y,Z)  from  the  NGS  program  gives  the  best  results  for  the 

larger  study  area.  For  practical  use  the  orthometric  heights  can  be  calculated  by 
utilizing  the  local  geoid  model. 

2.  The  areas  where  the  geoid  is  smoothest  lead  to  the  best  results.  This  was 
found  by  comparing  the  local  geoid  model  used  in  the  larger  Monterey  Bay  area  to 
the  local  geoid  model  used  on  the  NPS  campus.  Also  the  density  of  control  marks  on 
the  NPS  campus  was  much  higher  than  for  the  larger  study  area.  If  a  higher 
accuracy  is  required  in  a  large  area  where  the  geoid  may  be  undulating,  more 
control  marks  are  required.  These  control  marks  should  be  strategically  located 
throughout  the  network  in  order  to  determine  the  geoid  slope. 

3.  The  local  geoid  model  does  improve  the  accuracy  of  Rapp's  global  geoid 
model  locally.  The  accuracy  for  Rapp's  model  is  ±  2  m  in  an  area  55  km  x  55  km. 
The  accuracy  for  the  local  geoid  model  is  ±  2  cm  in  an  area  30  km  x  30  km  and  ±  1 
cm  in  an  area  0.5  km  x  0.5  km. 

4.  A  GPS  double  difference  fixed  solution  gives  the  best  ellipsoid  height 
differences,  Ah,  for  baseline  lengths  up  to  30  km. 

5.  The  standard  error  of  GPS  observations  of  the  ellipsoid  height  difference, 
Ah,  is  ±  2  cm  to  ±  4  cm,  so  there  is  no  need  to  perform  the  collocation  with  a  non¬ 
unit  weight  matrix.  Also  it  is  not  possible  to  reduce  the  error  of  the  GPS 
observation  below  ±  2  cm. 
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6.  There  is  no  need  for  local  surface  meteorological  corrections  to  the  GPS 
observations.  Default  meteorological  values  (1010  millibars,  20°  C  and  50% 

relative  humidity)  are  sufficient  to  meet  the  required  accuracy.  I 

7.  The  longer  the  duration  of  the  GPS  measurements,  the  better  the  result  will 
be.  The  accuracy  of  the  geoid  model  for  the  NPS  campus,  where  GPS 
measurements  were  made  over  four-  to  five-hour  periods,  is  better  than  the 

accuracy  of  the  geoid  height  model  in  the  larger  study  area,  where  measurements  I 

were  made  for  periods  of  only  ninety  minutes. 

8.  A  five-parameter  geoid  model  gives  the  best  solution  for  the  NPS  campus, 
while  a  six-parameter  geoid  model  gives  the  best  solution  for  the  larger  study  area. 

The  five-parameter  geoid  model  which  does  not  include  a  AZ  term  should  be  used  I 

with  caution  where  AZ  greatly  exceeds  5  m. 

Recommendations  are  as  follows: 

1 .  The  GPS  antenna  should  be  set  at  least  1  m  above  the  ground  to  reduce  the 

multiple  effects  from  the  ground.  The  height  of  the  antenna  should  be  carefully  | 

measured  before  and  after  the  data  collection. 

2.  The  configuration  of  the  control  marks  should  be  carefully  planned.  Sites 

should  be  selected  to  optimize  connections  to  stations  with  known  precise 
orthometric  height  differences,  minimizing  the  effects  of  obstructions.  Control  j 

marks  in  a  large  area  should  be  numerous  enough  to  meet  the  required  accuracy. 

3.  To  minimize  computation  errors  least  squares  adjustments  should  be  made 
using  double  precision  processing  when  solving  the  geoid  model.  For  large 

amounts  of  GPS  data  a  mainframe  computer  should  be  used  rather  than  a  j 

microcomputer  to  reduce  processing  time. 


I 


! 
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APPENDIX  A.  FORTRAN  PROGRAM  GPSCON 


********************************************************************** 


* 

* 

* 

* 

* 

* 

* 

* 


WRITTEN  BY  MA,  WEI-MING  05/08/88 

THIS  PROGRAM  READS  DATA  FROM  GEOID  DATA  FILE  AND  CALLS 
THE  SUBROUTINE  CONTOR  TO  PLOT  GEOID  HEIGHTS. 

THE  VARIABLES  USED  ARE: 


GEOHT 

NX 

NY 


GEOID  HEIGHTS  (m) 

NUMBER  OF  POINTS  IN  THE  X-DIRECTION 
NUMBER  OF  POINTS  IN  THE  Y-DIRECTION 


* 

* 

* 

* 

* 

* 

* 

* 


********************************************************************** 
REALM  GEOHT(22,15) 

DATA  NX,NY/22,15/ 


C 

C 

c 


READ  GEOID  HEIGHTS  FROM  GEOID  I  DATA  FILE 

CALL  EXCMSOFILEDEF  8  DISK  GEOID1  DATA  Al  ) 
DO  200  J  =  1,  NY 

DO  1001  =  NX,  1,-1 
READ(8,*)  GEOHT(I,J) 

100  CONTINUE 
200  CONTINUE 


C  PLOT  THE  GEOID  JGHT  BY  CALLING  SUBROUTINE  CONTOR 
C 

CALLCONT  ,R(GEOHT,NX,NY) 

STOP 

END 

SUBROUNTINE  CONTOR(A,NX,NY) 

********************************************************************** 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


WRITTEN  BY  MA,  WEI-MING  05/08/88 

THIS  SUBROUTINE  CONTOURS  AN  NX  BY  NY  ARRAY  OF 

REGULARLY  SPACED  POINTS. 

NOTE:  THIS  ARRAY  MUST  BE  REALM. 

THE  VARIABLES  USED  ARE: 

A  :  SINGLE  PRECISION  NX  BY  NY  ARRAY  OF 
REGULARLYSPACED  POINTS. 


NX 

NY 

ZINC 


NUMBER  OF  POINTS  IN  THE  X-DIRECTION 
NUMBER  OF  POINTS  IN  THE  Y-DIRECTION 
CONTOUR  INTERVAL 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


********************************************************************** 
DIMENSION  A(NX,NY) 

COMMON  WORK(SOOO) 

C 

C  SET  PARAMETERS  FOR  AXES 
C 
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CJUU  uuu  uuu  uuu  uu 


XORIG  =  -55.0 
XSTP  =  2.0 
XMAX  =  -34.0 
YORIG  =  35.0 
YSTP  =  2.0 
YMAX  =  49.0 


SET  CONTOUR  LEVEL 


ZINC  =  0.1 
CALL  COMPRS 
CALL  SETCLR(’CYAN’) 

SET  PAGE  AND  PLOT  SIZES,  SET  UP  AXES  FOR  PLOT 

CALL  PAGE(8.0,6.0) 

CALL  BCOMON(5000) 

CALL  AREA2D(7 .0,5.0) 

LABEL  AXES 

CALL  XNAME(’LONGITUDE  121  DEGREE  WEST  (MINUTES)$’,100) 
CALL  YNAME('LATTTUDE  36  DEGREE  NORTH  (MINUTES)$',100) 

CALL  GRAF(XORIG^STP^CMAX, YORIG, YSTP,  YMAX) 

CALL  FRAME 

MAKE  CONTOURS  AND  DRAW 

CALL  SETCLR('RED') 

CALL  CONMIN(3.0) 

CALL  CONANG(60.) 

CALL  CONLIN(0,’MYCON’, 'LABELS’,  1,10) 

CALL  CONMAK(A,NX,NY,ZINC) 

CALL  CONTURG, LABELS', 'DRAW) 

END  PLOT 

CALL  ENDPL(O) 

CALL  DONEPL 
RETURN 
END 

C  SUBROUTINE  MYCON(RARAY,LARAY) 

I********************************************************************** 
*  * 

*  THIS  SUBROUTINE  MAKES  NEGATIVE  CONTOURS  DASHED  AND  * 

*  THE  ZERO  LINE  HEAVYIER  * 

*  * 
*4^*******************************************************%*********** 

DIMENSION  RARA  Y  (2), LARA  Y  ( 1 ) 

CALL  RESET(’DASH') 

IF  (RARAY(l).GE.  0.)  GO  TO  10 
CALL  DASH 
10  RARAY(2)  =  1. 

IARAY(l)  =1 


IF  (RARAY(l).EQ.O.)  IARAY(1)  =  2 

RETURN 

END 
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APPENDIX  B.  BASIC  PROGRAM  LOBS.BASIC 


5  REM  LEAST  SQUARES  BY  OBSERVATION 

10  DIM  L(N%,1) 

20  DIM  A(N%,X%),  P(X%,X%),  EX(X%,1) 

30  DIM  CX(X%,1) 

40  DIM  AT(X%,N%),  R(50),  ATP(X%,X%) 

50  DIM  ATPA(X%,X%),  ATPL(X%,1),  AI(X%,X%) 

69  PRINT  "  INPUT  #  OF  OBSERVATIONS" 

70  NPUT  NO 

200  PRINT  "  INPUT  #  OF  PARAMETERS" 

270  INPUT  N1 

350  ERASE  L 

351  ERASE  A 

352  ERASE  AT 

353  ERASE  P 

354  ERASE  ATP 

355  ERASE  ATPL 

356  ERASE  AI 

357  ERASE  CX 

358  ERASE  EX 

359  ERASE  ATP  A 

380  DIM  L(NO,l),  A(NO,Nl),  AT(Nl,NO),  P(NO,NO),  ATP(Nl,NO), 

ATPL(N  1 , 1 ),  AI(N  1  ,N  1 ) 

390  DIM  CX(N  1,1),  EX  (NO,  1 ),  ATPA(N  1  ,N  1 ) 

400  ITERO  =  0 

1 1 80  PRINT  "  INPUT  COEFFICIENTS  OF  OBSERVATION  EQUATIONS  AND 
WEIGHTS" 

1190  FOR  K  =  1  TO  NO 
1200  FOR  J  =  1  TO  N1 
1220  PRINT  "COEFFICIENT  ",  K,  J 
1230  INPUT  A(K,J) 

1240  PRINT  A(K,J) 

1250  NEXT  J 

1260  PRINT  "  INPUT  OBSERVED  VALUE  ",  K 
1270  INPUT  L(K,1) 

1280  PRINT  "  INPUT  WEIGHT  OF  OBSERVED  VALUE  ",  K 
1300  INPUT  P(K,K) 

1390  NEXT  K 
1410  DOF  =  0 

1420  REM  LEAST  SQUARES 
1430  M  »  NO 
1450  L  =  N1 

1460  PRINT  "  M  =",  M,  "N1  =  ",  Nl,  "L  =  "L 
1470  FOR  I  =  1  TO  M 
1480  FOR  J  =  1  TO  L 
1490  AT(J,I)  =  A(I,J) 

1500  NEXT  J 
1510  NEXT  I 
1520  REM  ATP  =  AT  *  P 


1550  FOR  J  =1  TON 

1560  ATP(IJ)  =  0 

1570  FOR  K  =  1  TO  M 

1580  ATP(I,J)  =  ATP(I,J)  +  AT(U)  *  P(K,J) 

1590  NEXT  K 
1600  NEXT  J 
1610  NEXT  I 

1620  REM  ATPA  =  ATP  *  A 
1630  N  =  N1 
1640  FOR  1  =  1  TO  L 
1650  FOR  J  =  1  TO  N 
1660  ATPA(I,J)  =  0 
1680  ATPA(IJ)=0 
1690  FOR  K  =  1  TO  M 

17 10  ATPA(I,J)  =  ATPA(I,J)  +  ATP(I,K)  *  A(K,J) 
1720  NEXT  K 
1730  NEXT  J 
1740  NEXT  I 

1750  REM  ATPL  =  ATP  *  L 
1760  N  =  1 
1770  FOR  1=1  TO  L 
1780  FOR  J  =  1  TO  N 
1790  ATPL(IJ)  =  0 
1800  FOR  K  =  1  TOM 

1810  ATPL(IJ)  =  ATPL(I,J)  +  ATP(I,K)  *  L(K,J) 
1820  NEXT  K 
1830  NEXT  J 
1840  NEXT  I 

1970  REM  AI  =  INV(ATPA) 

1980  I  =  N1 

1990  M  =  N1 

2000  N  =  I  -  1 

2020  MI  =  M  -  1 

2030  FOR  J  =  1  TO  I 

2040  FOR  K  =  1  TO  I 

2050  AI(J,K)  =  ATPA(J,K) 

2055  NEXT  K 

2060  NEXT  J 

2070  FOR  K  =1  TO  I 

2080  FOR  J  =  1  TO  MI 

2090  R(J)  =  AI(1  J+l)  /  AI(  1,1) 

2100  NEXT  J 

2110  R(M)  =  1!  /  AI(1,1) 

2120  FOR  L  =  1  TO  N 
2130  FOR  J  =  1  TO  MI 

2140  AI(I,J)  =  AI(L+1,J+1)  -  AKL+1,1)  *  R(J) 
2150  NEXT  J 

2160  AI(L,M)  =  -AI(L+1,1)  *  R(M) 

2170  NEXT  L 
2180  FOR  J  =  1  TO  M 
2190  AI(I,J)  =  R(J) 

2200  NEXT  I 
2210  NEXT  K 
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2220  REM  CX  =  AI  *  ATPL 

2230  L  =  N1 

2240  N  =  1 

2250  M  =  N1 

2260  FOR  1=1  TO  L 

2270  FOR  J  =  1  TO  N 

2280  CX(I,J)  =  0 

2300  FOR  K  =  1  TO  M 

2310  CX(IJ)  =  CX(I,J)  +  AI(I,K)  *  ATPL  (K,J) 
2320  NEXT  K 

2325  PRINT  "  PARAMETER  #  ”,  I,  CX(I,J) 

2330  NEXT  J 

2340  NEXT  I 

2480  PRINT  "  RESIDUAL" 

2490  L  =  NO 

2500  N  =  1 

2510  M  =  N1 

2520  REM  EX  =  A  *  CX 

2530  FOR  1=1  TO  L 

2540  FOR  J  =  1  TO  N 

2560  EX(I,J)  =  0 

2570  FOR  K  =  1  TO  M 

2580  EX(IJ)  =  EX(I,J)  +  A(I,K)  *  CX(K,J) 

2590  NEXT  K 

2600  NEXT  J 

2610  NEXT  I 

2620  SD  =  0 

2630  FOR  K  =  1  TO  NO 

2640  V  *  EX(K,1)  -  L(K,1) 

2650  PRINT  K,  V 
2660  SD  =  SD  +  V  *  V 
2670  NEXT  K 

2680  SD  =  SQR(SD  /  (NO  -  N1  +  DOF)) 

2690  PRINT  SD 
2695  PRINT  "  STD.DEV",  SD 
2700  ITER  =  ITER  +  1 
2710  IF  ITER  <  3  GOTO  350 
2720  END 
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APPENDIX  C.  FORTRAN  PROGRAM  GPSDIS 


*********************************************************************** 


WRITTEN  BY  WEI-MING  MA  05/14/88 

THIS  PROGRAM  READS  DATA  FROM  TERMINAL  AND  COMPUTES 
THE  ORTHOMETRIC  HEIGHTS  FOR  THE  CHECK  MARKS  THEN 
COMPUTES  THE  DIFFERENCE. 

THE  VARIABLES  USED  ARE: 


RX,RY,RZ 
CX, CY.CZ 
DX.DY.DZ 
H 
LH 
DH 
GN 
LH 
DIF 

A.B.C.D 

A1.B1.C1 

D1.E1 

IPA.IPAl 

IW 


THE  COORDINATES  OF  THE  REFERENCE  MARK  (m) 
THE  COORDINATES  OF  THE  CHECK  MARK  (m) 

THE  COORDINATES  DIFFERENCE  (m) 

THE  ORTHOMETRIC  HEIGHT  OF  CHECK  MARK  (m) 
THE  LEVELED  ORTHOMETRIC  HEIGHT  (m) 

THE  ELLIPSOID  HEIGHT  DIFFERENCE  (m) 

THE  GLOBAL  GEOID  HEIGHT  OF  CHECK  MARK  (m) 
THE  LEVELED  ORTHOMETRIC  HEIGHT  (m) 

THE  DIFFERENCE  BETWEEN  H  AND  LH  (cm) 

THE  COEFFICIENTS  OF  5-PARAMETER  GEOID  MODEL 
THE  COEFFICIENTS  OF  6-PARAMETER  GEOID  MODEL 

THE  NUMBER  OF  PARAMETERS 
THE  OUTPUT  UNIT 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

*  * 
*********************************************************************** 

REAL* 8  RX,  RY,  RZ,  CX,  CY,  CZ,  DX,  DY,  DZ,  LH,  GN, 

REAL*8  A,  B,  C,  D,  Al,  Bl,  Cl,  Dl,  El,  H,  DH,  DIF, 

INTEGER  IPA.IPAl,  IW 

CHARACTER*  1  RESPN1,  RESPN2,  RESPN3 

IW  =  6 

IPA  =0 

IPA1  =0 


C 

C  READ  DATA  FROM  TERMINAL 
C 


10  PRINT  *,  'ENTER  #  OF  PARAMETERS' 

READ  *  IPA 

IF  (IPA  ,’eQ.  5  .OR.  IPA  .EQ.  6)  THEN 

PRINT  *,  'ENTER  THE  X  COORDINATE  IN  WGS  84  OF  THE  ', 
&  'REFERENCE  MARK’ 

READ  *,  RX 

PRINT  *,  'ENTER  THE  Y  COORDINATE  IN  WGS  84  OF  THE  ', 
&  'REFERENCE  MARK' 

READ  *,  RY 

PRINT  *,  'ENTER  THE  Z  COORDINATE  IN  WGS  84  OF  THE  ', 
&  'REFERENCE  MARK’ 

READ  *,  RZ 
ELSE 
STOP 
END  IF 
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IF  (RESPN1  .EQ.  'Y')  THEN 

PRINT  *,'DO  YOU  WANT  TO  CHANGE  THE  CHECH  MARK  ("Y"  OR  ', 

&  '"N")’ 

READ  *,  RESPN3 
IF  (RESPN3  .EQ.  ’Y')  GO  TO  30 
GO  TO  40 
END  IF 

20  PRINT  *,'DO  YOU  WANT  TO  CHANGE  THE  REFERENCE  MARKS  ("Y" 

&  ’OR  "N")' 

READ  *,  RESPNI 

IF  (RESPNI  .EQ.  'Y')  GO  TO  10 

30  PRINT  "VENTER  THE  X  COORDINATE  IN  WGS  84  OF  THE  CHECK.', 

&  ’MARK' 

READ  *  CX 

PRINT  VENTER  THE  Y  COORDINATE  IN  WGS  84  OF  THE  CHECK  MARK' 
READ  *  CY 

print  "■Venter  the  z  coordinate  in  wgs  84  of  the  check  mark' 

READ  *  CZ 

PRINT  "VENTER  THE  ELLIPSOID  DIFFERENCE,  DH' 

READ  *,  DH 

PRINT  "VENTER  THE  GLOBAL  GEOID  HEIGHT  OF  THE  CHECK  MARK', 

&  GN' 

READ  *,  GN 

PRINT  *, ENTER  THE  LEVELED  OTHOMETRIC  HEIGHT  OF  THE  CHECK  ’ 
&  ,  ’MARK’ 

READ  *,  LH 
DX  =  CX  -  RX 
DY  =  CY  -  RY 
DZ  =  CZ-  RZ 

IF  (IPA  .EQ.  IPA1  .AND.  RESPNI  .EQ.  ’N’  .AND.  IPA  .EQ.  5)  GO  TO  50 
IF  (IPA  .EQ.  IPA1  .AND.  RESPNI  .EQ.  ’N'  .AND.  IPA  .EQ.  6)  GO  TO  60 
40  IF  (IPA  .EQ.  5)  THEN 

PRINT  "V 'ENTER  THE  ELLIPSOID  HEIGHT  HO’ 

READ  *,  HO 

PRINT  *, 'ENTER  THE  COEFFICIENT  A' 

READ  *  A 

PRINT  *, 'ENTER  THE  COEFFICIENT  B' 

READ  *,  B 

PRINT  "VENTER  THE  COEFFICIENT  C' 

READ  *  C 

PRINT  *, 'ENTER  THE  COEFFICIENT  D' 

READ  *,  D 


THE  5-PARAMETER  GEOID  MODEL 


50  H  =  -GN  +  DH  +  HO  +  A*DY  +  B*DX"‘*2  +  C*DY**2  +  D*DX*DY 
ELSE 

PRINT  "VENTER  THE  ELLIPSOID  HEIGHT  HO' 

READ  *,  HO 

PRINT  "V ENTER  THE  COEFFICIENT  A’" 

READ  *  A1 

PRINT  *, 'ENTER  THE  COEFFICIENT  B’" 

READ  *,  B1 
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PRINT  *, 'ENTER  THE  COEFFICIENT  C"' 

READ  *,  Cl 

PRINT  *  ' 'ENTER  THE  COEFFICIENT  D'" 

READ  •  D1 

PRINT  *,'ENTER  THE  COEFFICIENT  E'" 

READ  *,  El 

THE  6-PARAMETER  GEOID  MODEL 

60  H  =  -GN  +  DH  +  HO  +  A1*DX  +  B1*DZ  +  C1*DX**2  +  D1*DY**2 
&  +  E1*DX*DY 

END  IF 

COMPUTE  THE  DIFFERENCE  AND  PRINT  THE  RESULT 

DIF  =  (  H  -  LH)  *  100. 

WRITE  (IW,  1)  IP  A 
WRITE(IW,  2)  LH,  H,  DIF 

1  FORMAT  (5X,  II, '  PARAMETERS') 

2  FORMAT  (/5X/H-LEVELING  =',  F10.3,  '  M' 

&  /5X,'H-GPS  =’,  F10.3,  '  M' 

&  /5X,'DIFFERENCE  =',  F10.3,  ’  CM'/) 

IPA1  =  IPA 
H  =  0. 

DIF  =  0. 

PRINT  *,  'MORE  COMPUTATION  ("Y"  OR  "N")' 

READ  *,  RESPN2 

IF  (RESPN2  .EQ.  'Y' )  GO  TO  20 

STOP 

END 
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APPENDIX  D.  FORTRAN  PROGRAM  DISTCO 


********************************************************************** 


* 


* 


WRITTEN  BY  MA,WEI_MING  05/11/88  * 

THIS  PROGRAM  COMPUTES  THE  DISTANCE  BETWEEN  THE  BENCH  * 
MARKS  THEN  CALLS  THE  SUBROUTINE  CTO  CALCULATE  THE  * 

COVARIANCE  BETWEEN  THE  BENCH  MARKS.  * 

THE  VARIABLES  USED  ARE:  * 

X,Y,Z  :  THE  COORDINATES  OF  X,Y,Z,  IN  WGS  84  (m)  * 

DIS  :  THE  DISTANCES  BETWEEN  THE  MARKS  (m)  * 

CQ  :  THE  COVARIANCES  BETWEEN  THE  BENCH  MARKS  * 


********************************************************************** 
REAL*8  X(10),  Y(10),  Z(10) 

REAL*8  DIS(10,10),  CQ(10,10) 

DATA  DIS/100*0./,  CQ/100*0./ 

rw  =  9 

N  =  7 
C 

C  READ  X,  Y,  Z  FROM  THE  DATA  FILE 
C 

CALL  EXCMSOFILEDEF  8  DISK  MBXYZ  DATA  Al') 

DO  100  I  =  1,  N 

READ(8,*)  X(I),  Y(I),  Z(I)) 

100  CONTINUE 
C 

C  COMPUTE  THE  DISTANCE  AND  PRINT  THE  RESULTS 
C 

DO  300  I  =  1,  N 
WRITE(IW,1)  I 

1  FORMAT(/ IX, 'FROM  BENCH  MARK  #',  12/  TO',  'DISTANCE  (m)’/ 

&  IX,  53('-')) 

DO  200  J  =  1,  N 

IF  (J  .EQ.  I)  GO  TO  10 

DIS(IJ)  =  SQRT((X(J)  -  X(I))  **  2  +  (Y(J)  -  Y(I))  **  2 
&  +(Z(J)-Z(I))**2) 

io  WRrrE(iw,2)  j,  dis(U) 

2  FORMAT!  11X,  12,  14X,  G20.8) 

200  CONTINUE 

300  CONTINUE 
C 

C  CALL  SUBROUTINE  C  TO  CALCULATE  THE  CQ  AND  PRINT  THE  RESULTS 
C 

WRITE  (TW, 3) 

3  FORMAT(/3X,  'CQ  =  ') 

CALL  QDIS.N.CQ) 

DO  400  J  =  1 ,  N 

WRITE0W.4)  (CQ(I,J),  I  =  1,  N) 

4  FORMAT(3X,  7(G16.7,  2X) 

400  CONTINUE 


STOP 

END 

SUBROUNTINE  C(DIS,N,CQ) 

*  *  *  *  *  *  *  *  *  ****  *  *  *  *  Hi*  *  *  *  *  %  *  *  *  1)1*  *  *  *  *  **  *  *  *  4c  **  *  *  *  *  **  *  *  *  *  **  *  **  *  *  *  *  *  *  *  *  *  *  *  *  *  * 


WRITTEN  BY  MA,WEI_MING  05/11/88  * 

THIS  SUBROUTINE  COMPUTES  THE  COVARIANCES  OF  THE  RANDOM  * 
ERROR  DISTRIBUTION  FUNCTION:  * 

C(R)  =  A  +  B  *  SIN(D*DIS)  * 

THE  VARIABLES  USED  ARE:  * 

A  :  STANDARD  DEVIATION  OF  THE  OBSERVATIONS  * 

B  :  CONSTANT  * 

C  :  CONSTANT  (RADIANS/METERS)  * 

DIS  :  DISTANCE  (METERS)  * 

* 

id*  ***  *  *****  *  *  *  *  ****  *  *****  *  *****  *  **  **  *  *  ****  *  ***  ****  4c  *  *************  *  ***  * 

REAL*8  DIS(10,10),  CQ(10,10),  A,  B,  D 
A  =  0.137153 
B  =  0.147951 
D  =  2.85958 
DO  200  J  =  1,  N 
DO  100  I  =  1,  N 

CQ(IJ)  =  A  -  B  *  SIN(D  *  DIS(IJ)  /  33300.) 

100  CONTINUE 
200  CONTINUE 
RETURN 
END 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 
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